---
title: "Replication Code"
subtitle: "Taking the Cloth: Social norms and elite cues increase support for masks among white Evangelical Americans"
author: "Claire Adida, Christina Cottiero, Leonardo Falabella, Isabel Gotti, ShahBano Ijaz, Gregoire Phillips, Michael Seese"
date: "Last Updated: `r format(Sys.time(), '%d %B %Y')`"
output: 
  html_document:
      toc: true
      toc_depth: 4
      toc_float:
        collapsed: false
      code_folding: hide
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, warning = FALSE, message = FALSE)

library("tidyverse")
library("modelsummary")
library("cobalt")
library("AER")
library("mediation")
library("MASS")
library("ggeffects")
library("emmeans")
library("esc")
library("factoextra")
```

# Overview

The following document contains replication code for the analyses included in the manuscript and supplementary information of the paper *Taking the Cloth: Social norms and elite cues increase support for masks among white Evangelical Americans*.

## Data & Code
All data generated by this study are contained in `working.rds`, available through the Journal of Experimental Political Science Repository at the [Harvard Dataverse](https://dataverse.harvard.edu/dataverse/xps).  See the Supplementary Information document for a list and description of all included variables.

Annotated replication code is provided in this document; R scripts are also available through Harvard Dataverse.  See the included `README.txt` file and  Supplementary Information (SI) document for additional notes and information. 

```{r, data}
## Loading the data
working <- readRDS(file = "working.rds")

## Recode religiosity index to correct a data transfer error from Lucid
working$rel1x <- dplyr::recode(working$rel1,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel2x <- dplyr::recode(working$rel2,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel3x <- dplyr::recode(working$rel3,
                               `1` = 1,
                               `2` = 6,
                               `3` = 2,
                               `4` = 3,
                               `5` = 4,
                               `6` = 5)

working$rel.ix <- working$rel1x + working$rel2x + working$rel3x
```

# Manuscript

This section provides replication code for our main results, as well as descriptive statistics we reference in the manuscript. Replication code for results in the manuscript for which we point the reader to the Supplementary Information are included in the [following section](#si).

## Descriptive Statistics

In the **Results** section of the paper (Section 3), we provide descriptive statistics on the gender, age, and partisan makeup of our sample.  We provide the full set of descriptive statistics in our SI document (replication code below).

### Gender
```{r gender_descriptives, fig.width = 10}
ggplot(data = working, aes(x = factor(gender))) + 
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5, size = 2.5) +
  labs(title = "Gender Distribution of Sample",
       x = "Gender",
       y = "Frequency") +
  theme_bw()
```

### Age
```{R age_descriptives, fig.width = 10}
ggplot(data = working, aes(x = age)) + 
  geom_histogram(aes(y = ..density..)) +
  geom_vline(aes(xintercept = mean(age, na.rm = TRUE)),
             color ="red", linetype = "dashed", size = 0.5) +
  labs(title = "Age Distribution of Sample",
       x = "Age",
       y = "Density",
       caption = "Mean Age = 54.23") +
  theme_bw()
```

### Party ID
```{R party_descriptives, fig.width = 10}
ggplot(data = working, aes(x = factor(party.3))) + 
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5, size = 2.5) +
  labs(title = "Partisan Distribution of Sample", 
       subtitle = "Three Categories",
       x = "Party ID",
       y = "Frequency") +
  theme_bw()
```

## Average Treatment Effects{#m-ate}

We test the two main hypotheses using the pre-registered specification:
$$y_i = \alpha + \beta{T_i} + \gamma{X_i} + \epsilon_i$$
Where $y_i$ is the outcome of interest for respondent $i$, $\beta$ is the estimated treatment effect, $T_i$ is the assignment to treatment, and $X_i$ is a vector of controls.  The following code estimates this model and reproduces **Figure 1** in the manuscript.

```{R main_results, fig.width = 10}
## Subset the working dataset to include only respondents eligible to click petition link
w2 <- working %>% 
  filter(behave > 4)

## Specify model with all controls
x <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"

## Run the models
c1 <- lm(as.formula(paste("therm", x)), data = working)
c2 <- lm(as.formula(paste("p1.wash", x)), data = working)
c3 <- lm(as.formula(paste("p2.dist", x)), data = working)
c4 <- lm(as.formula(paste("p3.mask", x)), data = working)
c5 <- lm(as.formula(paste("p4.pray", x)), data = working)
c6 <- lm(as.formula(paste("behave", x)), data = working)
c7 <- lm(as.formula(paste("clicked", x)), data = w2)

## Create DF for coefficient plot
x1 <- tidy(c1, conf.int = TRUE)
x2 <- tidy(c2, conf.int = TRUE)
x3 <- tidy(c3, conf.int = TRUE)
x4 <- tidy(c4, conf.int = TRUE)
x5 <- tidy(c5, conf.int = TRUE)
x6 <- tidy(c6, conf.int = TRUE)
x7 <- tidy(c7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(c1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(c2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(c3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(c4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(c5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(c6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(c7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
ggplot(y) +
  geom_hline(yintercept = 0, lty = "11", colour = "grey30") +
  geom_errorbar(aes(cond, ymin = conf.low, ymax = conf.high), width = 0.2) +
  geom_errorbar(aes(cond, ymin = low, ymax = high), lwd = 1.15, width = 0) +
  geom_point(aes(cond, estimate)) +
  labs(colour = "", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out)) + 
  labs(title = "Coefficient Plot: OLS Model with Controls",
       caption = "Shown with 90% and 95% confidence intervals.
       Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

```

We contextualize these results by estimating the average treatment effect in terms of standard deviations, which we calculate using the `esc` package.

```{R ate_sd}
## Norms
esc_B(b = 0.45198, sdy = 2.995311, grp1n = 597, grp2n = 592, es.type = "g")

## Endorsement
esc_B(b = 0.340044, sdy = 2.995311, grp1n = 597, grp2n = 591, es.type = "g")
```

## Prior Knowledge

In the manuscript, we cite statistics indicating that respondents were less likely to have prior knowledge of the information provided by our norms treatment than the information provided by our endorsement treatment.  Below, we show that only 49.15% of respondents in the norms condition knew the information provided by the vignette, compared to 81.38% of respondents in the endorsement condition.

```{R knowledge_of_treatment, fig.width = 10}
working2 <- working %>% 
  unite("knew", 13:15, na.rm = TRUE, remove = FALSE)

working2$knew <- as.numeric(working2$knew)

working2$knewit <- dplyr::recode(working2$knew,
                                 `3` = "Knew It",
                                 `2` = "Didn't Know It",
                                 `1` = "Don't Believe It",
                                 .default = NA_character_)

working2$knewit <- as.factor(working2$knewit)

t <- table(working2$treat, working2$knewit)
t
prop.table(t, 1)

ggplot(data = working2, aes(x = factor(knewit))) +
  geom_bar(stat = "count") +
  geom_text(stat = "count", aes(label = ..count..), vjust = -0.5, size = 2.5) +
  labs(title = "Prior Knowledge of Treatment Information, by Treatment Condition",
       x = "Prior Knowledge",
       y = "Frequency") +
  facet_wrap(~ treat) +
  theme_bw()
```


## Principal Component Analysis

In Footnote 4, we state: "A principal components analysis on our six non-behavioral outcome measures shows that a single dimension explains a relatively high proportion of the overall variation in our outcome measures—about 78.28%.'' The following code replicates this PCA, and generates a scree plot and a variable contribution plot.  Because the semi-behavioral outcome ("would click") is not closely related to the other outcome variables, we drop this variable from the synthetic index described [later in this documet](#syndex) 

```{R pca, fig.width = 10}
# PCA For Footnote 4 ###########################################################
## Simplify the dataset
x <- working %>% dplyr::select(
  treat, therm, p1.wash, p2.dist, p3.mask, p4.pray, behave)

x2 <- x %>% dplyr::select(
  "Thermometer" = therm, 
  "Hand Washing" = p1.wash, 
  "Distancing" = p2.dist, 
  "Mask Wearing" = p3.mask, 
  "Prayer" = p4.pray, 
  "Would Click" = behave)

## Run the PCA
pc1 <- prcomp(na.omit(x2), scale = FALSE)
summary(pc1)


## PCA Plots ###################################################################
fviz_eig(pc1)

fviz_pca_var(pc1, 
                     # col.var = "contrib", 
                     repel = TRUE)
```


## Other Results

The manuscript alludes to the results of a compliance analysis and an exploratory analysis that looks at partisan identification, both of which are contained in the SI.  Code to replicate these results can be found in the [subsequent section](#si) of this document.

# Supplementrary Information{#si}

This section contains replication code for all tables, figures, and analyses contained in the SI document.

## Descriptive Statistics

### Covariates

```{R ds_covariates_1}
x1 <- working %>% dplyr::select(
  `Age` = age,
  `Follows News` = news.follow,
  `Follows Covid` = covid.follow,
  `Covid Threat` = covid.threat)

datasummary_skim(x1)
```

### Covariates (Categorical)

```{R ds_covariates_2}
x2 <- working %>% dplyr::select(
  `Gender` = gender,
  `Education` = edu,
  `Party ID` = party.id,
  `News Source` = news.source,
  `Covid News Source` = covid.source,
  `Tested for Covid` = covid.tested,
  `Friends with Covid` = covid.friends)

datasummary_skim(x2, type = "categorical")
```

### Religious Variables

```{R ds_religious_vars}
x3 <- working %>% dplyr::select(
  `Regularly Attends Church` = rel1x,
  `Spiritual Values` = rel2x,
  `Americans Should be More Religious` = rel3x,
  `Religiosity Index` = rel.ix,
  `Religiosity (Self-Reported)` = religiosity,
  `Religious Motivation - Insurance` = c1, 
  `Religious Motivation - Morality` = c2, 
  `Religious Motivation - Social Capital` = c3, 
  `Religious Motivation - Heuristics` = c4,
  `Religious Motivation - Tradition` = c5, 
  `Religious Motivation - Social Identity` = c6, 
  `Religion Offers Comfort` = ei1,
  `Takes Religious Approach to Life` = ei2,
  `Enjoys Social Aspects of Church` = ei3,
  `Extrinsic / Intrinsic Motivation Index` = ei.index)

datasummary_skim(x3)

rm(x1, x2, x3, summaryname)
```


## Balance

### Covariate Balance

```{R bal_covariates}
dem <- working %>% dplyr::select(
  `Condition` = treat,
  `Age` = age,
  `Gender` = gender,
  `Education` = edu,
  `Party ID` = party.id,
  `Follows News` = news.follow,
  `News Source` = news.source,
  `Follows Covid` = covid.follow,
  `Covid News Source` = covid.source,
  `Covid Threat` = covid.threat,
  `Tested for Covid` = covid.tested,
  `Friends with Covid` = covid.friends)

datasummary_balance(~Condition,
                    data = dem,
                    fmt = '%.2f')
```

### Religious Variables Balance

```{R bal_religious_vars}
mech <- working %>% dplyr::select(
  `Condition` = treat,
  `Regularly Attends Church` = rel1x,
  `Spiritual Values` = rel2x,
  `Americans Should be More Religious` = rel3x,
  `Religiosity Index` = rel.ix,
  `Religiosity (Self-Reported)` = religiosity,
  `Religious Motivation - Insurance` = c1, 
  `Religious Motivation - Morality` = c2, 
  `Religious Motivation - Social Capital` = c3, 
  `Religious Motivation - Heuristics` = c4,
  `Religious Motivation - Tradition` = c5, 
  `Religious Motivation - Social Identity` = c6, 
  `Religion Offers Comfort` = ei1,
  `Takes Religious Approach to Life` = ei2,
  `Enjoys Social Aspects of Church` = ei3,
  `Extrinsic / Intrinsic Motivation Index` = ei.index)

datasummary_balance(~Condition,
                    data = mech,
                    fmt = '%.2f')
```

### Love Plot

```{R love_plot, fig.width = 10, fig.height=8}
## Data Cleaning
dem <- dem[complete.cases(dem), ]
dem2 <- subset(dem, select = -c(Condition))

bt <- bal.tab(dem2, treat = dem$Condition, m.threshold = .1, 
              multi.summary = TRUE)

## Clean up variable names in Excel and save as .csv
# var.names(bt, type = "df", minimal = FALSE)

nn <- read.csv("loveplot_newnames.csv")

## Love Plot
love.plot(dem2, treat = dem$Condition, 
                thresholds = c(m = .1), 
                binary = "std",
                which.treat = .all, 
                var.names = nn,
                abs = FALSE) + 
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  labs(subtitle = "Standardized Mean Differences Across Comparison Groups") +
  theme(legend.position = "none")

## Cleaning Up
rm(list=setdiff(ls(), "working"))
```


## Correlates of Mask Wearing
In the manuscript, before we present our main analyses, we provide a brief summary of the respondents in our sample who were most likely to wear a mask.  This analysis considers only respondents in the control condition; the dependent variable is our thermometer score.  We find that older, more educated, and less Republican respondents are more likely to agree with the importance of mask-wearing in public.

### Demographic Correlates

```{R corr_masking_dem_1}
control <- working %>%
  subset(treat == "Control")

control$ag <- factor(control$age.group, ordered = FALSE)
control$eg <- factor(control$edu, ordered = FALSE)


### OLS Models #################################################################
## OLS Model: Demographics
x1 <- lm(therm ~ gender + ag + eg + party.f + news.follow + covid.follow + 
           covid.threat + covid.tested + covid.friends, data = control)
# summary(x1)

## OLS Table
cm <- c("(Intercept)" = "Intercept",
        "genderMale" = "Male",
        "ag30-39" = "Age 30-39",
        "ag40-49" = "Age 40-49",
        "ag50-59" = "Age 50-59",
        "ag60-69" = "Age 60-69",
        "ag70+" = "Age 70+",
        "egHigh School" = "High School",
        "egSome College" = "Some College",
        "egVocational School" = "Vocational School",
        "egAssociate's" = "Associate's",
        "egBachelor's" = "Bachelor's",
        "egSome Graduate School" = "Some Graduate School",
        "egMaster's" = "Master's",
        "egProfessional Degree" = "Professional Degree",
        "egDoctorate" = "Doctorate",
        "party.fDemocrat" = "Democrat",
        "party.fLean Democrat" = "Lean Democrat",
        "party.fIndependent" = "Independent",
        "party.fLean Republican" = "Lean Republican",
        "party.fRepublican" = "Republican",
        "party.fStrong Republican" = "Strong Republican",
        "news.follor" = "Follows News",
        "covid.threat" = "Covid Threat",
        "covid.testedYes" = "Tested for Covid",
        "covid.friendsYes" = "Friends with Covid")

modelsummary(x1,
             # output = "latex",
             title = "Demographic Correlates of Mask Wearing",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')
```

The following charts plot the predicted thermometer score at each value of age group, education, and partisanship.  Note that these predictions are based on the above model, which contains the full set of demographic covariates

```{R corr_masking_dem_2, fig.width = 10}
## Linear Predictions by Demographic
ggplot(ggpredict(x1, terms = "ag"), aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Age Group", y = "Predicted Value")

ggplot(ggpredict(x1, terms = "eg"), aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Education", y = "Predicted Value")

ggplot(ggpredict(x1, terms = "party.f"), aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Party", y = "Predicted Value")

rm("x1", "cm")
```

### Religious Correlates

Below, we estimate a set of models that look at the religious correlates of mask-wearing.

```{R corr_masking_rel, fig.width = 10}
## Religious Correlates of Mask Wearing ########################################
### Data Cleaning ##############################################################
## Extrinsic / Intrinsic Motivations
control$ei1.f <- factor(control$ei1)
control$ei2.f <- factor(control$ei2)
control$ei3.f <- factor(control$ei3)

control$ei1.f <- dplyr::recode(control$ei1.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

control$ei2.f <- dplyr::recode(control$ei2.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

control$ei3.f <- dplyr::recode(control$ei3.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree")

## Religious Motivations
control$c.all <- NA
control$c.all[which(control$c1 == 1)] <- 1
control$c.all[which(control$c2 == 1)] <- 2
control$c.all[which(control$c3 == 1)] <- 3
control$c.all[which(control$c4 == 1)] <- 4
control$c.all[which(control$c5 == 1)] <- 5
control$c.all[which(control$c6 == 1)] <- 6

# table(control$c.all)
control$c.all <- factor(control$c.all)

## Religiosity
control$religiosity.f <- factor(control$religiosity)
control$religiosity.f <- dplyr::recode(control$religiosity.f,
                                       "1" = "Anti-Religious",
                                       "2" = "Not at all Religious",
                                       "3" = "Slightly Religious",
                                       "4" = "Moderately Religious",
                                       "5" = "Very Religious")

## Religiosity Index
control$rel1.f <- factor(control$rel1)
control$rel2.f <- factor(control$rel2)
control$rel3.f <- factor(control$rel3)

control$rel1.f <- dplyr::recode(control$rel1.f,
                                "1" = "Strongly Disagree",
                                "3" = "Disagree",
                                "4" = "Somewhat Disagree",
                                "5" = "Somewhat Agree",
                                "6" = "Agree",
                                "2" = "Strongly Agree")

control$rel2.f <- dplyr::recode(control$rel2.f,
                                "1" = "Strongly Disagree",
                                "3" = "Disagree",
                                "4" = "Somewhat Disagree",
                                "5" = "Somewhat Agree",
                                "6" = "Agree",
                                "2" = "Strongly Agree")

control$rel3.f <- dplyr::recode(control$rel3.f,
                                "1" = "Strongly Disagree",
                                "3" = "Disagree",
                                "4" = "Somewhat Disagree",
                                "5" = "Somewhat Agree",
                                "6" = "Agree",
                                "2" = "Strongly Agree")


### OLS Models #################################################################
## OLS Model: Extrinsic / Intrinsic Motivations
x1 <- lm(therm ~ ei1.f + ei2.f + ei3.f, data = control)
# summary(x1)

## OLS Model: Religious Motivations
x2 <- lm(therm ~ c.all, data = control)
# summary(x2)

## OLS Model: Religiosity
x3 <- lm(therm ~ religiosity.f, data = control)
# summary(x3)

## OLS Model: Religiosity Index
x4 <- lm(therm ~ rel1.f + rel2.f + rel3.f, data = control)
# summary(x4)

## OLS Model: All Religious Variables
x5 <- lm(therm ~ ei1.f + ei2.f + ei3.f + c.all + religiosity.f + rel1.f + 
           rel2.f + rel3.f, data = control)
# summary(x5)

## OLS Table
cm <- c("(Intercept)" = "Intercept",
        "ei1.fDisagree" = "Religion Offers Comfort: Disagree",
        "ei1.fNeither Agree nor Disagree" = "Religion Offers Comfort: Neither Agree nor Disagree",
        "ei1.fAgree" = "Religion Offers Comfort: Agree",
        "ei1.fStrongly Agree" = "Religion Offers Comfort: Strongly Agree",
        "ei2.fDisagree" = "Religious Approach to Life: Disagree",
        "ei2.fNeither Agree nor Disagree" = "Religious Approach to Life: Neither Agree nor Disagree",
        "ei2.fAgree" = "Religious Approach to Life: Agree",
        "ei2.fStrongly Agree" = "Religious Approach to Life: Strongly Agree",
        "ei3.fDisagree" = "Enjoys Social Aspects of Church: Disagree",
        "ei3.fNeither Agree nor Disagree" = "Enjoys Social Aspects of Church: Neither Agree nor Disagree",
        "ei3.fAgree" = "Enjoys Social Aspects of Church: Agree",
        "ei3.fStrongly Agree" = "Enjoys Social Aspects of Church: Strongly Agree",
        "c.all2" = "Morality",
        "c.all3" = "Social Capital",
        "c.all4" = "Heuristics",
        "c.all5" = "Tradition",
        "c.all6" = "Social Identity",
        "religiosity.fSlightly Religious" = "Religiosity (SR) Slightly Religious",
        "religiosity.fModerately Religious" = "Religiosity (SR) Moderately Religious",
        "religiosity.fVery Religious" = "Religiosity (SR) Very Religious",
        "rel1.fStrongly Agree" = "Regularly Attends Church: Strongly Agree",
        "rel1.fDisagree" = "Regularly Attends Church: Disagree",
        "rel1.fSomewhat Disagree" = "Regularly Attends Church: Somewhat Disagree",
        "rel1.fSomewhat Agree" = "Regularly Attends Church: Somewhat Agree",
        "rel1.fAgree" = "Regularly Attends Church: Agree",
        "rel2.fStrongly Agree" = "Spiritual Values: Strongly Agree",
        "rel2.fDisagree" = "Spiritual Values: Disagree",
        "rel2.fSomewhat Disagree" = "Spiritual Values: Somewhat Disagree",
        "rel2.fSomewhat Agree" = "Spiritual Values: Somewhat Agree",
        "rel2.fAgree" = "Spiritual Values: Agree",
        "rel3.fStrongly Agree" = "Americans Should be More Religious: Strongly Agree",
        "rel3.fDisagree" = "Americans Should be More Religious: Disagree",
        "rel3.fSomewhat Disagree" = "Americans Should be More Religious: Somewhat Disagree",
        "rel3.fSomewhat Agree" = "Americans Should be More Religious: Somewhat Agree",
        "rel3.fAgree" = "Americans Should be More Religious: Agree")

modelsummary(list(x1, x2, x3, x4, x5),
             # output = "latex",
             title = "Religious Correlates of Mask Wearing",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

rm("x1", "x2", "x3", "x4", "x5", "cm")
```

### Partisan Correlates

We also regress respondents' thermometer scores on party ID to assess the naive correlation between masking-wearing and partisanship.  The predicted thermometer scores by party ID are shown in the plots below.

```{R corr_masking_party, fig.width = 10}
## Party ID and Mask Wearing ###################################################
### Data Cleaning ##############################################################
control$party.f2 <- factor(control$party.f, levels = c("Strong Democrat", 
                                                       "Democrat", 
                                                       "Lean Democrat",
                                                       "Independent", 
                                                       "Lean Republican", 
                                                       "Republican", 
                                                       "Strong Republican"))

control$party.32 <- factor(control$party.3, levels = c("Democrat",
                                                       "Independent", 
                                                       "Republican"))


### OLS Models #################################################################
## 7 Category Party ID
x1 <- lm(therm ~ party.f2, data = control)
# summary(x1)

m.x1 <- ggpredict(x1, terms = "party.f2")

ggplot(m.x1, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.1) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Party ID", y = "Predicted Value")

## 3 Category Party ID
x2 <- lm(therm ~ party.32, data = control)
# summary(x2)

m.x2 <- ggpredict(x2, terms = "party.32")

ggplot(m.x2, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin = conf.low, ymax = conf.high), width = 0.075) +
  geom_text(aes(label = round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 2) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Linear Predictions",
       subtitle = "Thermometer (Control Condition)",
       x  = "Party ID", y = "Predicted Value")

## Cleaning Up
rm(x1, m.x1, x2, m.x2, control)
```

## Average Treatment Effects

In this section, we re-estimate the ATE models for each outcome variable (see [above](#m-ate)), and generate the OLS tables that appear in the SI.  For aesthetic reasons and ease of interpretability, we omit controls from the following tables, though the full set of results can be recovered from the replication code.

```{R ate}
# Average Treatment Effects: OLS Tables ########################################
## Data Cleaning ###############################################################
w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"


## OLS Models ##################################################################
## Baseline Models
a1 <- lm(therm ~ treat, data = working)
a2 <- lm(p1.wash ~ treat, data = working)
a3 <- lm(p2.dist ~ treat, data = working)
a4 <- lm(p3.mask ~ treat, data = working)
a5 <- lm(p4.pray ~ treat, data = working)
a6 <- lm(behave ~ treat, data = working)
a7 <- lm(clicked ~ treat, data = w2)

## Demographic Controls
b1 <- lm(as.formula(paste("therm", iv1)), data = working)
b2 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
b3 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
b4 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
b5 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
b6 <- lm(as.formula(paste("behave", iv1)), data = working)
b7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

## All Controls
c1 <- lm(as.formula(paste("therm", iv2)), data = working)
c2 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
c3 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
c4 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
c5 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
c6 <- lm(as.formula(paste("behave", iv2)), data = working)
c7 <- lm(as.formula(paste("clicked", iv2)), data = w2)


## OLS Tables ##################################################################
## Group models by outcome variable
therm <- list(a1, b1, c1)
wash <- list(a2, b2, c2)
dist <- list(a3, b3, c3)
mask <- list(a4, b4, c4)
pray <- list(a5, b5, c5)
beh <- list(a6, b6, c6)
click <- list(a7, b7, c7)

## Coefficients
cm = c("(Intercept)" = "Intercept",
       "treatEndorsement" = "Endorsement",
       "treatNorms" = "Norms")

## Adding Rows to the Tables
row <- tribble(~x0, ~x1, ~x2, ~x3,
               "Controls", 
               "No", "Demographic", "All" )
attr(row, 'position') <- c(7)


## Tables by Outcome Variable
modelsummary(therm,
             title = "Thermometer",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(wash,
             title = "Hand Washing",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(dist,
             title = "Distancing",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(mask,
             title = "Mask Wearing",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(pray,
             title = "Prayer",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(beh,
             title = "Would Click",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

modelsummary(click,
             title = "Clicked",
             coef_map = cm,
             add_rows = row,
             stars = TRUE,
             gof_omit = 'IC|Log|Adj')

## Cleaning Up
rm(list=setdiff(ls(), "working"))
```

## Compliance Average Causal Effects

For our pre-registered compliance analysis, we asked each respondent whether or not they already knew the information provided by the vignette, and whether or not they believed it.  Those who responded "I do not believe this" were coded as non-compliers. 

The following table shows the breakdown of compliant vs. non-compliant respondents in each treatment condition.

```{R cace_1}
# Compliance Average Causal Effects ############################################
## Generate Compliance Variables & Subset Data #################################
## Compliance Variables
# as.data.frame(colnames(working))

working2 <- working %>% 
  unite("knew", 13:15, na.rm = TRUE, remove = FALSE)

working2$knew <- as.numeric(working2$knew)

working2$knewit <- dplyr::recode(working2$knew,
                                 `3` = "Knew It",
                                 `2` = "Didn't Know It",
                                 `1` = "Don't Believe It",
                                 .default = NA_character_)

working2$knewit <- as.factor(working2$knewit)
# table(working2$knewit)

working2$comply <- dplyr::recode(working2$knewit,
                                 "Don't Believe It" = "Non-Compliant",
                                 "Didn't Know It" = "Compliant",
                                 "Knew It" = "Compliant")

# table(working2$treat)
table(working2$treat, working2$comply)
```

Using a pre-registered instrumental variable approach, we estimate the treatment effect for compliant respondents for each outcome variable.

```{R cace_2}
## Endorsement Dataset
e <- subset(working2, treat != "Norms")
e$treat.e <- factor(e$treat, 
                    levels = c("Control", "Endorsement"),
                    ordered = TRUE)

e$t.e[e$treat.e == "Control"] <- 0
e$t.e[e$treat.e == "Endorsement"] <- 1
e$tu.e[e$comply == "Non-Compliant"] <- 0
e$tu.e[e$comply == "Compliant"] <- 1
e$tu2.e <- e$tu.e
e$tu2.e[e$t.e == 0] <- 0
# Coding Control as all non-compliant
# Because you can't comply with a treatment you didn't get
# tu2.e is now the compliance variable: = 1 if compliant | endorsement condition

e2 <- e %>% 
  filter(behave > 4)

## Norms Dataset
n <- subset(working2, treat != "Endorsement")
n$treat.n <- factor(n$treat, 
                    levels = c("Control", "Norms"),
                    ordered = TRUE)

n$t.n[n$treat.n == "Control"] <- 0
n$t.n[n$treat.n == "Norms"] <- 1
n$tu.n[n$comply == "Non-Compliant"] <- 0
n$tu.n[n$comply == "Compliant"] <- 1
n$tu2.n <- n$tu.n
n$tu2.n[n$t.n == 0] <- 0
# tu2.n is now the compliance variabe: = 1 if compliant | norms condition

n2 <- n %>% 
  filter(behave > 4)


## Endorsement IV Models #######################################################
## Baseline
e1.1 <- ivreg(therm ~ tu2.e, ~ t.e, data = e)
e1.2 <- ivreg(p1.wash ~ tu2.e, ~ t.e, data = e)
e1.3 <- ivreg(p2.dist ~ tu2.e, ~ t.e, data = e)
e1.4 <- ivreg(p3.mask ~ tu2.e, ~ t.e, data = e)
e1.5 <- ivreg(p4.pray ~ tu2.e, ~ t.e, data = e)
e1.6 <- ivreg(behave ~ tu2.e, ~ t.e, data = e)
e1.7 <- ivreg(clicked ~ tu2.e, ~ t.e, data = e2)


## Demographic Controls
e2.1 <- ivreg(therm ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.2 <- ivreg(p1.wash ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.3 <- ivreg(p2.dist ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.4 <- ivreg(p3.mask ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.5 <- ivreg(p4.pray ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.6 <- ivreg(behave ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e)
e2.7 <- ivreg(clicked ~ tu2.e, ~ t.e + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = e2)

## All Controls
e3.1 <- ivreg(therm ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.1 <- ivreg(therm ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.2 <- ivreg(p1.wash ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.3 <- ivreg(p2.dist ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.4 <- ivreg(p3.mask ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.5 <- ivreg(p4.pray ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.6 <- ivreg(behave ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e)
e3.7 <- ivreg(clicked ~ tu2.e, ~ t.e + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = e2)


## Norms IV Models #############################################################
## Baseline
n1.1 <- ivreg(therm ~ tu2.n, ~ t.n, data = n)
n1.2 <- ivreg(p1.wash ~ tu2.n, ~ t.n, data = n)
n1.3 <- ivreg(p2.dist ~ tu2.n, ~ t.n, data = n)
n1.4 <- ivreg(p3.mask ~ tu2.n, ~ t.n, data = n)
n1.5 <- ivreg(p4.pray ~ tu2.n, ~ t.n, data = n)
n1.6 <- ivreg(behave ~ tu2.n, ~ t.n, data = n)
n1.7 <- ivreg(clicked ~ tu2.n, ~ t.n, data = n2)

## Demographic Controls
n2.1 <- ivreg(therm ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.2 <- ivreg(p1.wash ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.3 <- ivreg(p2.dist ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.4 <- ivreg(p3.mask ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.5 <- ivreg(p4.pray ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.6 <- ivreg(behave ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n)
n2.7 <- ivreg(clicked ~ tu2.n, ~ t.n + gender + age + education + party.f + covid.threat + covid.tested + covid.friends, data = n2)

## All Controls
n3.1 <- ivreg(therm ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.1 <- ivreg(therm ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.2 <- ivreg(p1.wash ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.3 <- ivreg(p2.dist ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.4 <- ivreg(p3.mask ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.5 <- ivreg(p4.pray ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.6 <- ivreg(behave ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n)
n3.7 <- ivreg(clicked ~ tu2.n, ~ t.n + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends, data = n2)


## CACE Tables #################################################################
therm <- list(e1.1, e2.1, e3.1, n1.1, n2.1, n3.1)
wash <- list(e1.2, e2.2, e3.2, n1.2, n2.2, n3.2)
dist <- list(e1.3, e2.3, e3.3, n1.3, n2.3, n3.3)
mask <- list(e1.4, e2.4, e3.4, n1.4, n2.4, n3.4)
pray <- list(e1.5, e2.5, e3.5, n1.5, n2.5, n3.5)
beh <- list(e1.6, e2.6, e3.6, n1.6, n2.6, n3.6)
click <- list(e1.7, e2.7, e3.7, n1.7, n2.7, n3.7)

row <- tribble(~x0, ~x1, ~x2, ~x3, ~x4, ~x5, ~x6,
               "Controls", 
               "No", "Demographic", "All", 
               "No", "Demographic", "All" )
attr(row, 'position') <- c(7)

cm = c("(Intercept)" = "Intercept",
       "tu2.e" = "Endorsement (Compliant)",
       "tu2.n" = "Norms (Compliant)")

modelsummary(therm,
             title = "Thermometer",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(wash,
             title = "Hand Washing",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(dist,
             title = "Distancing",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(mask,
             title = "Mask Wearing",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(pray,
             title = "Prayer",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(beh,
             title = "Would Click",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

modelsummary(click,
             title = "Clicked",
             coef_map = cm, 
             add_rows = row,
             stars = TRUE)

## Cleaning Up
rm(list=setdiff(ls(), "working"))
```

## Mediation Analyses

Sections 2.4.2 and 2.4.3 in our pre-analysis plan register a set of mediation analyses to test a religiosity mechanism and religious motivation mechanism.  Consistent with our PAP, we estimate the following models using the full set of controls.  The following plots show the average causal mediation effects (ACME), the average direct effects (ADE), and the combined (total) effects.  Note that the following code chunks may take some time to run.

### Religiosity (Self Reported)

```{R med_religiosity, fig.width = 10, fig.height = 10, cache = TRUE}
# Mediation Analyses ###########################################################
x <- " ~ treat + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Religiosity Mechanism #######################################################
working <- rename(working, "religiosity.raw" = "religiosity")

working$religiosity <- factor(working$religiosity.raw)
working$religiosity <- dplyr::recode(working$religiosity,
                                     "1" = "Anti-Religious",
                                     "2" = "Not at all Religious",
                                     "3" = "Slightly Religious",
                                     "4" = "Moderately Religious",
                                     "5" = "Very Religious")

working2 <- working %>% 
  filter(behave > 4)

## Mediator
h4.med <- polr(as.formula(paste("religiosity", x)), data = working, Hess = TRUE)
h4.med2 <- polr(as.formula(paste("religiosity", x)), data = working2, Hess = TRUE)

## Outcomes
h4.therm <- lm(as.formula(paste("therm", x, " + religiosity")), data = working)
h4.p1 <- lm(as.formula(paste("p1.wash", x, " + religiosity")), data = working)
h4.p2 <- lm(as.formula(paste("p2.dist", x, " + religiosity")), data = working)
h4.p3 <- lm(as.formula(paste("p3.mask", x, " + religiosity")), data = working)
h4.p4 <- lm(as.formula(paste("p4.pray", x, " + religiosity")), data = working)
h4.behave <- lm(as.formula(paste("behave", x, " + religiosity")), data = working)
h4.click <- lm(as.formula(paste("clicked", x, " + religiosity")), data = working2)

## Endorsement
h4.e.therm <- mediate(h4.med, h4.therm, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

h4.e.p1 <- mediate(h4.med, h4.p1, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

h4.e.p2 <- mediate(h4.med, h4.p2, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

h4.e.p3 <- mediate(h4.med, h4.p3, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

h4.e.p4 <- mediate(h4.med, h4.p4, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

h4.e.behave <- mediate(h4.med, h4.behave, treat = "treat", mediator = "religiosity", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

h4.e.click <- mediate(h4.med2, h4.click, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## Norms
h4.n.therm <- mediate(h4.med, h4.therm, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

h4.n.p1 <- mediate(h4.med, h4.p1, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 1000)

h4.n.p2 <- mediate(h4.med, h4.p2, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 1000)

h4.n.p3 <- mediate(h4.med, h4.p3, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 1000)

h4.n.p4 <- mediate(h4.med, h4.p4, treat = "treat", mediator = "religiosity", 
                   dropobs = TRUE, control.value = "Control", 
                   treat.value  = "Norms", robustSE = TRUE, sims = 1000)

h4.n.behave <- mediate(h4.med, h4.behave, treat = "treat", mediator = "religiosity", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

h4.n.click <- mediate(h4.med2, h4.click, treat = "treat", mediator = "religiosity", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## Mediation Plots
## Endorsement
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(h4.e.therm, main = "Thermometer")
plot(h4.e.p1, main = "Policy: Hand Washing")
plot(h4.e.p2, main = "Policy: Distancing")
plot(h4.e.p3, main = "Policy: Mask Wearing")
plot(h4.e.p4, main = "Policy: Prayer")
plot(h4.e.behave, main = "Would Click")
plot(h4.e.click, main = "Clicked")
mtext("Mediator Variable: Religiosity (Self-Reported) \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Norms
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(h4.n.therm, main = "Thermometer")
plot(h4.n.p1, main = "Policy: Hand Washing")
plot(h4.n.p2, main = "Policy: Distancing")
plot(h4.n.p3, main = "Policy: Mask Wearing")
plot(h4.n.p4, main = "Policy: Prayer")
plot(h4.n.behave, main = "Would Click")
plot(h4.n.click, main = "Clicked")
mtext("Mediator Variable: Religiosity (Self-Reported) \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)
```

### Religious Motivations

```{R med_rel_motiv, fig.width = 10, fig.height = 10, cache = TRUE}
## Religious Motivations #######################################################
rm(list=setdiff(ls(), c("working", "working2", "x")))

## Mediator
mc1 <- lm(paste("c1", x), data = working)
mc2 <- lm(paste("c2", x), data = working)
mc3 <- lm(paste("c3", x), data = working)
mc4 <- lm(paste("c4", x), data = working)
mc5 <- lm(paste("c5", x), data = working)
mc6 <- lm(paste("c6", x), data = working)

mc1.2 <- lm(paste("c1", x), data = working2)
mc2.2 <- lm(paste("c2", x), data = working2)
mc3.2 <- lm(paste("c3", x), data = working2)
mc4.2 <- lm(paste("c4", x), data = working2)
mc5.2 <- lm(paste("c5", x), data = working2)
mc6.2 <- lm(paste("c6", x), data = working2)

## Outcomes
## Motivation 1
o1c1 <- lm(paste("therm", x, " + c1"), data = working)
o2c1 <- lm(paste("p1.wash", x, " + c1"), data = working)
o3c1 <- lm(paste("p2.dist", x, " + c1"), data = working)
o4c1 <- lm(paste("p3.mask", x, " + c1"), data = working)
o5c1 <- lm(paste("p4.pray", x, " + c1"), data = working)
o6c1 <- lm(paste("behave", x, " + c1"), data = working)
o7c1 <- lm(paste("clicked", x, " + c1"), data = working2)

## Motivation 2
o1c2 <- lm(paste("therm", x, " + c2"), data = working)
o2c2 <- lm(paste("p1.wash", x, " + c2"), data = working)
o3c2 <- lm(paste("p2.dist", x, " + c2"), data = working)
o4c2 <- lm(paste("p3.mask", x, " + c2"), data = working)
o5c2 <- lm(paste("p4.pray", x, " + c2"), data = working)
o6c2 <- lm(paste("behave", x, " + c2"), data = working)
o7c2 <- lm(paste("clicked", x, " + c2"), data = working2)

## Motivation 3
o1c3 <- lm(paste("therm", x, " + c3"), data = working)
o2c3 <- lm(paste("p1.wash", x, " + c3"), data = working)
o3c3 <- lm(paste("p2.dist", x, " + c3"), data = working)
o4c3 <- lm(paste("p3.mask", x, " + c3"), data = working)
o5c3 <- lm(paste("p4.pray", x, " + c3"), data = working)
o6c3 <- lm(paste("behave", x, " + c3"), data = working)
o7c3 <- lm(paste("clicked", x, " + c3"), data = working2)

## Motivation 4
o1c4 <- lm(paste("therm", x, " + c4"), data = working)
o2c4 <- lm(paste("p1.wash", x, " + c4"), data = working)
o3c4 <- lm(paste("p2.dist", x, " + c4"), data = working)
o4c4 <- lm(paste("p3.mask", x, " + c4"), data = working)
o5c4 <- lm(paste("p4.pray", x, " + c4"), data = working)
o6c4 <- lm(paste("behave", x, " + c4"), data = working)
o7c4 <- lm(paste("clicked", x, " + c4"), data = working2)

## Motivation 5
o1c5 <- lm(paste("therm", x, " + c5"), data = working)
o2c5 <- lm(paste("p1.wash", x, " + c5"), data = working)
o3c5 <- lm(paste("p2.dist", x, " + c5"), data = working)
o4c5 <- lm(paste("p3.mask", x, " + c5"), data = working)
o5c5 <- lm(paste("p4.pray", x, " + c5"), data = working)
o6c5 <- lm(paste("behave", x, " + c5"), data = working)
o7c5 <- lm(paste("clicked", x, " + c5"), data = working2)

## Motivation 6
o1c6 <- lm(paste("therm", x, " + c6"), data = working)
o2c6 <- lm(paste("p1.wash", x, " + c6"), data = working)
o3c6 <- lm(paste("p2.dist", x, " + c6"), data = working)
o4c6 <- lm(paste("p3.mask", x, " + c6"), data = working)
o5c6 <- lm(paste("p4.pray", x, " + c6"), data = working)
o6c6 <- lm(paste("behave", x, " + c6"), data = working)
o7c6 <- lm(paste("clicked", x, " + c6"), data = working2)

## Mediation
## Motivation 1, Endorsement Condition
m1o1c1 <- mediate(mc1, o1c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o2c1 <- mediate(mc1, o2c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o3c1 <- mediate(mc1, o3c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o4c1 <- mediate(mc1, o4c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o5c1 <- mediate(mc1, o5c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o6c1 <- mediate(mc1, o6c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o7c1 <- mediate(mc1.2, o7c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

## Motivation 1, Norms Condition
m2o1c1 <- mediate(mc1, o1c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o2c1 <- mediate(mc1, o2c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o3c1 <- mediate(mc1, o3c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o4c1 <- mediate(mc1, o4c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o5c1 <- mediate(mc1, o5c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o6c1 <- mediate(mc1, o6c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o7c1 <- mediate(mc1.2, o7c1, treat = "treat", mediator = "c1", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

## Motivation 2, Endorsement Condition
m1o1c2 <- mediate(mc2, o1c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o2c2 <- mediate(mc2, o2c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o3c2 <- mediate(mc2, o3c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o4c2 <- mediate(mc2, o4c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o5c2 <- mediate(mc2, o5c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o6c2 <- mediate(mc2, o6c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o7c2 <- mediate(mc2.2, o7c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

## Motivation 2, Norms Condition
m2o1c2 <- mediate(mc2, o1c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o2c2 <- mediate(mc2, o2c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o3c2 <- mediate(mc2, o3c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o4c2 <- mediate(mc2, o4c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o5c2 <- mediate(mc2, o5c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o6c2 <- mediate(mc2, o6c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o7c2 <- mediate(mc2.2, o7c2, treat = "treat", mediator = "c2", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

## Motivation 3, Endorsement Condition
m1o1c3 <- mediate(mc3, o1c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o2c3 <- mediate(mc3, o2c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o3c3 <- mediate(mc3, o3c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o4c3 <- mediate(mc3, o4c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o5c3 <- mediate(mc3, o5c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o6c3 <- mediate(mc3, o6c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o7c3 <- mediate(mc3.2, o7c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

## Motivation 3, Norms Condition
m2o1c3 <- mediate(mc3, o1c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o2c3 <- mediate(mc3, o2c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o3c3 <- mediate(mc3, o3c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o4c3 <- mediate(mc3, o4c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o5c3 <- mediate(mc3, o5c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o6c3 <- mediate(mc3, o6c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o7c3 <- mediate(mc3.2, o7c3, treat = "treat", mediator = "c3", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

## Motivation 4, Endorsement Condition
m1o1c4 <- mediate(mc4, o1c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o2c4 <- mediate(mc4, o2c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o3c4 <- mediate(mc4, o3c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o4c4 <- mediate(mc4, o4c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o5c4 <- mediate(mc4, o5c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o6c4 <- mediate(mc4, o6c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o674 <- mediate(mc4.2, o7c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

## Motivation 4, Norms Condition
m2o1c4 <- mediate(mc4, o1c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o2c4 <- mediate(mc4, o2c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o3c4 <- mediate(mc4, o3c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o4c4 <- mediate(mc4, o4c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o5c4 <- mediate(mc4, o5c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o6c4 <- mediate(mc4, o6c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o7c4 <- mediate(mc4.2, o7c4, treat = "treat", mediator = "c4", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

## Motivation 5, Endorsement Condition
m1o1c5 <- mediate(mc5, o1c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o2c5 <- mediate(mc5, o2c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o3c5 <- mediate(mc5, o3c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o4c5 <- mediate(mc5, o4c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o5c5 <- mediate(mc5, o5c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o6c5 <- mediate(mc5, o6c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o7c5 <- mediate(mc5.2, o7c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

## Motivation 5, Norms Condition
m2o1c5 <- mediate(mc5, o1c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o2c5 <- mediate(mc5, o2c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o3c5 <- mediate(mc5, o3c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o4c5 <- mediate(mc5, o4c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o5c5 <- mediate(mc5, o5c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o6c5 <- mediate(mc5, o6c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o7c5 <- mediate(mc5.2, o7c5, treat = "treat", mediator = "c5", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

## Motivation 6, Endorsement Condition
m1o1c6 <- mediate(mc6, o1c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o2c6 <- mediate(mc6, o2c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o3c6 <- mediate(mc6, o3c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o4c6 <- mediate(mc6, o4c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o5c6 <- mediate(mc6, o5c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o6c6 <- mediate(mc6, o6c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

m1o7c6 <- mediate(mc6.2, o7c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Endorsement" , robustSE = TRUE, sims = 1000)

## Motivation 6, Norms Condition
m2o1c6 <- mediate(mc6, o1c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o2c6 <- mediate(mc6, o2c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o3c6 <- mediate(mc6, o3c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o4c6 <- mediate(mc6, o4c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o5c6 <- mediate(mc6, o5c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o6c6 <- mediate(mc6, o6c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

m2o7c6 <- mediate(mc6.2, o7c6, treat = "treat", mediator = "c6", dropobs = TRUE, 
                  control.value = "Control", treat.value  = "Norms" , robustSE = TRUE, sims = 1000)

## Mediation Plots
## Motivation 1, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c1, main = "Thermometer")
plot(m1o2c1, main = "Policy: Hand Washing")
plot(m1o3c1, main = "Policy: Distancing")
plot(m1o4c1, main = "Policy: Mask Wearing")
plot(m1o5c1, main = "Policy: Prayer")
plot(m1o6c1, main = "Would Click")
plot(m1o7c1, main = "Clicked")
mtext("Mediator Variable: Insurance \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 1, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c1, main = "Thermometer")
plot(m2o2c1, main = "Policy: Hand Washing")
plot(m2o3c1, main = "Policy: Distancing")
plot(m2o4c1, main = "Policy: Mask Wearing")
plot(m2o5c1, main = "Policy: Prayer")
plot(m2o6c1, main = "Would Click")
plot(m2o7c1, main = "Clicked")
mtext("Mediator Variable: Insurance \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 2, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c2, main = "Thermometer")
plot(m1o2c2, main = "Policy: Hand Washing")
plot(m1o3c2, main = "Policy: Distancing")
plot(m1o4c2, main = "Policy: Mask Wearing")
plot(m1o5c2, main = "Policy: Prayer")
plot(m1o6c2, main = "Would Click")
plot(m1o7c2, main = "Clicked")
mtext("Mediator Variable: Morality \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 2, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c2, main = "Thermometer")
plot(m2o2c2, main = "Policy: Hand Washing")
plot(m2o3c2, main = "Policy: Distancing")
plot(m2o4c2, main = "Policy: Mask Wearing")
plot(m2o5c2, main = "Policy: Prayer")
plot(m2o6c2, main = "Would Click")
plot(m2o7c2, main = "Clicked")
mtext("Mediator Variable: Morality \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 3, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c3, main = "Thermometer")
plot(m1o2c3, main = "Policy: Hand Washing")
plot(m1o3c3, main = "Policy: Distancing")
plot(m1o4c3, main = "Policy: Mask Wearing")
plot(m1o5c3, main = "Policy: Prayer")
plot(m1o6c3, main = "Would Click")
plot(m1o7c3, main = "Clicked")
mtext("Mediator Variable: Social Capital \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 3, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c3, main = "Thermometer")
plot(m2o2c3, main = "Policy: Hand Washing")
plot(m2o3c3, main = "Policy: Distancing")
plot(m2o4c3, main = "Policy: Mask Wearing")
plot(m2o5c3, main = "Policy: Prayer")
plot(m2o6c3, main = "Would Click")
plot(m2o7c3, main = "Clicked")
mtext("Mediator Variable: Social Capital \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 4, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c4, main = "Thermometer")
plot(m1o2c4, main = "Policy: Hand Washing")
plot(m1o3c4, main = "Policy: Distancing")
plot(m1o4c4, main = "Policy: Mask Wearing")
plot(m1o5c4, main = "Policy: Prayer")
plot(m1o6c4, main = "Would Click")
plot(m1o674, main = "Clicked")
mtext("Mediator Variable: Heuristics \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 4, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c4, main = "Thermometer")
plot(m2o2c4, main = "Policy: Hand Washing")
plot(m2o3c4, main = "Policy: Distancing")
plot(m2o4c4, main = "Policy: Mask Wearing")
plot(m2o5c4, main = "Policy: Prayer")
plot(m2o6c4, main = "Would Click")
plot(m2o7c4, main = "Clicked")
mtext("Mediator Variable:  Heuristics \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 5, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c5, main = "Thermometer")
plot(m1o2c5, main = "Policy: Hand Washing")
plot(m1o3c5, main = "Policy: Distancing")
plot(m1o4c5, main = "Policy: Mask Wearing")
plot(m1o5c5, main = "Policy: Prayer")
plot(m1o6c5, main = "Would Click")
plot(m1o7c5, main = "Clicked")
mtext("Mediator Variable: Tradition \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 5, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c5, main = "Thermometer")
plot(m2o2c5, main = "Policy: Hand Washing")
plot(m2o3c5, main = "Policy: Distancing")
plot(m2o4c5, main = "Policy: Mask Wearing")
plot(m2o5c5, main = "Policy: Prayer")
plot(m2o6c5, main = "Would Click")
plot(m2o7c5, main = "Clicked")
mtext("Mediator Variable: Tradition \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 6, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m1o1c6, main = "Thermometer")
plot(m1o2c6, main = "Policy: Hand Washing")
plot(m1o3c6, main = "Policy: Distancing")
plot(m1o4c6, main = "Policy: Mask Wearing")
plot(m1o5c6, main = "Policy: Prayer")
plot(m1o6c6, main = "Would Click")
plot(m1o7c6, main = "Clicked")
mtext("Mediator Variable: Social Identity \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Motivation 6, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(m2o1c6, main = "Thermometer")
plot(m2o2c6, main = "Policy: Hand Washing")
plot(m2o3c6, main = "Policy: Distancing")
plot(m2o4c6, main = "Policy: Mask Wearing")
plot(m2o5c6, main = "Policy: Prayer")
plot(m2o6c6, main = "Would Click")
plot(m2o7c6, main = "Clicked")
mtext("Mediator Variable: Social Identity \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Cleaning Up
rm(list=setdiff(ls(), c("working", "x")))
```

### Extrinsic / Intrinsic Motivations

```{R med_ei_motiv, fig.width = 10, fig.height = 10, cache = TRUE}
## Extrinsic / Intrinsic Motivations ###########################################
## Data Cleaning
working$ei1.f <- factor(working$ei1)
working$ei2.f <- factor(working$ei2)
working$ei3.f <- factor(working$ei3)

working$ei1.f <- dplyr::recode(working$ei1.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

working$ei2.f <- dplyr::recode(working$ei2.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

working$ei3.f <- dplyr::recode(working$ei3.f,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Neither Agree nor Disagree",
                               "4" = "Agree",
                               "5" = "Strongly Agree") 

working2 <- working %>% 
  filter(behave > 4)

## Mediator
m.ei1 <- polr(as.formula(paste("ei1.f", x)), data = working, Hess = TRUE)
m.ei2 <- polr(as.formula(paste("ei2.f", x)), data = working, Hess = TRUE)
m.ei3 <- polr(as.formula(paste("ei3.f", x)), data = working, Hess = TRUE)

m2.ei1 <- polr(as.formula(paste("ei1.f", x)), data = working2, Hess = TRUE)
m2.ei2 <- polr(as.formula(paste("ei2.f", x)), data = working2, Hess = TRUE)
m2.ei3 <- polr(as.formula(paste("ei3.f", x)), data = working2, Hess = TRUE)

## Outcomes
## EI 1 - Religion Offers Comfort
ei1.therm <- lm(as.formula(paste("therm", x, " + ei1.f")), data = working)
ei1.p1 <- lm(as.formula(paste("p1.wash", x, " + ei1.f")), data = working)
ei1.p2 <- lm(as.formula(paste("p2.dist", x, " + ei1.f")), data = working)
ei1.p3 <- lm(as.formula(paste("p3.mask", x, " + ei1.f")), data = working)
ei1.p4 <- lm(as.formula(paste("p4.pray", x, " + ei1.f")), data = working)
ei1.behave <- lm(as.formula(paste("behave", x, " + ei1.f")), data = working)
ei1.click <- lm(as.formula(paste("clicked", x, " + ei1.f")), data = working2)

## EI 2 - Religious Approach to Life
ei2.therm <- lm(as.formula(paste("therm", x, " + ei2.f")), data = working)
ei2.p1 <- lm(as.formula(paste("p1.wash", x, " + ei2.f")), data = working)
ei2.p2 <- lm(as.formula(paste("p2.dist", x, " + ei2.f")), data = working)
ei2.p3 <- lm(as.formula(paste("p3.mask", x, " + ei2.f")), data = working)
ei2.p4 <- lm(as.formula(paste("p4.pray", x, " + ei2.f")), data = working)
ei2.behave <- lm(as.formula(paste("behave", x, " + ei2.f")), data = working)
ei2.click <- lm(as.formula(paste("clicked", x, " + ei2.f")), data = working2)

## EI 3 - Enjoys Social Aspects of Church
ei3.therm <- lm(as.formula(paste("therm", x, " + ei3.f")), data = working)
ei3.p1 <- lm(as.formula(paste("p1.wash", x, " + ei3.f")), data = working)
ei3.p2 <- lm(as.formula(paste("p2.dist", x, " + ei3.f")), data = working)
ei3.p3 <- lm(as.formula(paste("p3.mask", x, " + ei3.f")), data = working)
ei3.p4 <- lm(as.formula(paste("p4.pray", x, " + ei3.f")), data = working)
ei3.behave <- lm(as.formula(paste("behave", x, " + ei3.f")), data = working)
ei3.click <- lm(as.formula(paste("clicked", x, " + ei3.f")), data = working2)

## Mediation
## EI 1,  Endorsement Condition
ei1.e.therm <- mediate(m.ei1, ei1.therm, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei1.e.p1 <- mediate(m.ei1, ei1.p1, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei1.e.p2 <- mediate(m.ei1, ei1.p2, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei1.e.p3 <- mediate(m.ei1, ei1.p3, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei1.e.p4 <- mediate(m.ei1, ei1.p4, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei1.e.behave <- mediate(m.ei1, ei1.behave, treat = "treat", mediator = "ei1.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei1.e.click <- mediate(m2.ei1, ei1.click, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## EI 1, Norms Condition
ei1.n.therm <- mediate(m.ei1, ei1.therm, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei1.n.p1 <- mediate(m.ei1, ei1.p1, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei1.n.p2 <- mediate(m.ei1, ei1.p2, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei1.n.p3 <- mediate(m.ei1, ei1.p3, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei1.n.p4 <- mediate(m.ei1, ei1.p4, treat = "treat", mediator = "ei1.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei1.n.behave <- mediate(m.ei1, ei1.behave, treat = "treat", mediator = "ei1.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei1.n.click <- mediate(m2.ei1, ei1.click, treat = "treat", mediator = "ei1.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## EI2, Endorsement Condition
ei2.e.therm <- mediate(m.ei2, ei2.therm, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei2.e.p1 <- mediate(m.ei2, ei2.p1, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei2.e.p2 <- mediate(m.ei2, ei2.p2, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei2.e.p3 <- mediate(m.ei2, ei2.p3, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei2.e.p4 <- mediate(m.ei2, ei2.p4, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei2.e.behave <- mediate(m.ei2, ei2.behave, treat = "treat", mediator = "ei2.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei2.e.click <- mediate(m2.ei2, ei2.click, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## EI 2, Norms Condition
ei2.n.therm <- mediate(m.ei2, ei2.therm, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei2.n.p1 <- mediate(m.ei2, ei2.p1, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei2.n.p2 <- mediate(m.ei2, ei2.p2, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei2.n.p3 <- mediate(m.ei2, ei2.p3, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei2.n.p4 <- mediate(m.ei2, ei2.p4, treat = "treat", mediator = "ei2.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei2.n.behave <- mediate(m.ei2, ei2.behave, treat = "treat", mediator = "ei2.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei2.n.click <- mediate(m2.ei2, ei2.click, treat = "treat", mediator = "ei2.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## EI 3, Endorsement Condition
ei3.e.therm <- mediate(m.ei3, ei3.therm, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei3.e.p1 <- mediate(m.ei3, ei3.p1, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei3.e.p2 <- mediate(m.ei3, ei3.p2, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei3.e.p3 <- mediate(m.ei3, ei3.p3, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei3.e.p4 <- mediate(m.ei3, ei3.p4, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei3.e.behave <- mediate(m.ei3, ei3.behave, treat = "treat", mediator = "ei3.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

ei3.e.click <- mediate(m2.ei3, ei3.click, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## EI 3, Norms Condition
ei3.n.therm <- mediate(m.ei3, ei3.therm, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei3.n.p1 <- mediate(m.ei3, ei3.p1, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei3.n.p2 <- mediate(m.ei3, ei3.p2, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei3.n.p3 <- mediate(m.ei3, ei3.p3, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei3.n.p4 <- mediate(m.ei3, ei3.p4, treat = "treat", mediator = "ei3.f", 
                    dropobs = TRUE, control.value = "Control", 
                    treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei3.n.behave <- mediate(m.ei3, ei3.behave, treat = "treat", mediator = "ei3.f", 
                        dropobs = TRUE, control.value = "Control", 
                        treat.value  = "Norms", robustSE = TRUE, sims = 1000)

ei3.n.click <- mediate(m2.ei3, ei3.click, treat = "treat", mediator = "ei3.f", 
                       dropobs = TRUE, control.value = "Control", 
                       treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## Mediation Plots
## EI 1, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei1.e.therm, main = "Thermometer")
plot(ei1.e.p1, main = "Policy: Hand Washing")
plot(ei1.e.p2, main = "Policy: Distancing")
plot(ei1.e.p3, main = "Policy: Mask Wearing")
plot(ei1.e.p4, main = "Policy: Prayer")
plot(ei1.e.behave, main = "Would Click")
plot(ei1.e.click, main = "Clicked")
mtext("Mediator Variable: Religion Offers Comfort \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 1, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei1.n.therm, main = "Thermometer")
plot(ei1.n.p1, main = "Policy: Hand Washing")
plot(ei1.n.p2, main = "Policy: Distancing")
plot(ei1.n.p3, main = "Policy: Mask Wearing")
plot(ei1.n.p4, main = "Policy: Prayer")
plot(ei1.n.behave, main = "Would Click")
plot(ei1.n.click, main = "Clicked")
mtext("Mediator Variable: Religion Offers Comfort \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 2, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei2.e.therm, main = "Thermometer")
plot(ei2.e.p1, main = "Policy: Hand Washing")
plot(ei2.e.p2, main = "Policy: Distancing")
plot(ei2.e.p3, main = "Policy: Mask Wearing")
plot(ei2.e.p4, main = "Policy: Prayer")
plot(ei2.e.behave, main = "Would Click")
plot(ei2.e.click, main = "Clicked")
mtext("Mediator Variable: Religious Approach to Life \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 2, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei2.n.therm, main = "Thermometer")
plot(ei2.n.p1, main = "Policy: Hand Washing")
plot(ei2.n.p2, main = "Policy: Distancing")
plot(ei2.n.p3, main = "Policy: Mask Wearing")
plot(ei2.n.p4, main = "Policy: Prayer")
plot(ei2.n.behave, main = "Would Click")
plot(ei2.n.click, main = "Clicked")
mtext("Mediator Variable: Religious Approach to Life \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 3, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei3.e.therm, main = "Thermometer")
plot(ei3.e.p1, main = "Policy: Hand Washing")
plot(ei3.e.p2, main = "Policy: Distancing")
plot(ei3.e.p3, main = "Policy: Mask Wearing")
plot(ei3.e.p4, main = "Policy: Prayer")
plot(ei3.e.behave, main = "Would Click")
plot(ei3.e.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Endorsement Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## EI 3, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(ei3.n.therm, main = "Thermometer")
plot(ei3.n.p1, main = "Policy: Hand Washing")
plot(ei3.n.p2, main = "Policy: Distancing")
plot(ei3.n.p3, main = "Policy: Mask Wearing")
plot(ei3.n.p4, main = "Policy: Prayer")
plot(ei3.n.behave, main = "Would Click")
plot(ei3.n.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Norms Condition", 
      outer=TRUE,  side = 3, cex=1, line=-0.5)

## Cleaning Up
rm(list=setdiff(ls(), c("working", "x")))
```

### Religiosity Index

```{R med_relindex_motiv, fig.width = 10, fig.height = 10, cache = TRUE}
## Religiosity Index ###########################################################
working$rel1x <- factor(working$rel1x)
working$rel2x <- factor(working$rel2x)
working$rel3x <- factor(working$rel3x)

working$rel1x <- dplyr::recode(working$rel1x,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Somewhat Disagree",
                               "4" = "Somewhat Agree",
                               "5" = "Agree",
                               "6" = "Strongly Agree")

working$rel2x <- dplyr::recode(working$rel2x,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Somewhat Disagree",
                               "4" = "Somewhat Agree",
                               "5" = "Agree",
                               "6" = "Strongly Agree")

working$rel3x <- dplyr::recode(working$rel3x,
                               "1" = "Strongly Disagree",
                               "2" = "Disagree",
                               "3" = "Somewhat Disagree",
                               "4" = "Somewhat Agree",
                               "5" = "Agree",
                               "6" = "Strongly Agree")

working2 <- working %>% 
  filter(behave > 4)

## Mediator
m.rel1x <- polr(as.formula(paste("rel1x", x)), data = working, Hess = TRUE)
m.rel2x <- polr(as.formula(paste("rel2x", x)), data = working, Hess = TRUE)
m.rel3x <- polr(as.formula(paste("rel3x", x)), data = working, Hess = TRUE)

m2.rel1x <- polr(as.formula(paste("rel1x", x)), data = working2, Hess = TRUE)
m2.rel2x <- polr(as.formula(paste("rel2x", x)), data = working2, Hess = TRUE)
m2.rel3x <- polr(as.formula(paste("rel3x", x)), data = working2, Hess = TRUE)

## Outcomes
## Religiosity Index 1 - Regularly Attends Church
rel1x.therm <- lm(as.formula(paste("therm", x, " + rel1x")), data = working)
rel1x.p1 <- lm(as.formula(paste("p1.wash", x, " + rel1x")), data = working)
rel1x.p2 <- lm(as.formula(paste("p2.dist", x, " + rel1x")), data = working)
rel1x.p3 <- lm(as.formula(paste("p3.mask", x, " + rel1x")), data = working)
rel1x.p4 <- lm(as.formula(paste("p4.pray", x, " + rel1x")), data = working)
rel1x.behave <- lm(as.formula(paste("behave", x, " + rel1x")), data = working)
rel1x.click <- lm(as.formula(paste("clicked", x, " + rel1x")), data = working2)

## Religiosity Index 2 - Spiritual Values
rel2x.therm <- lm(as.formula(paste("therm", x, " + rel2x")), data = working)
rel2x.p1 <- lm(as.formula(paste("p1.wash", x, " + rel2x")), data = working)
rel2x.p2 <- lm(as.formula(paste("p2.dist", x, " + rel2x")), data = working)
rel2x.p3 <- lm(as.formula(paste("p3.mask", x, " + rel2x")), data = working)
rel2x.p4 <- lm(as.formula(paste("p4.pray", x, " + rel2x")), data = working)
rel2x.behave <- lm(as.formula(paste("behave", x, " + rel2x")), data = working)
rel2x.click <- lm(as.formula(paste("clicked", x, " + rel2x")), data = working2)

## Religiosity Index 3 - Enjoys Social Aspects
rel3x.therm <- lm(as.formula(paste("therm", x, " + rel3x")), data = working)
rel3x.p1 <- lm(as.formula(paste("p1.wash", x, " + rel3x")), data = working)
rel3x.p2 <- lm(as.formula(paste("p2.dist", x, " + rel3x")), data = working)
rel3x.p3 <- lm(as.formula(paste("p3.mask", x, " + rel3x")), data = working)
rel3x.p4 <- lm(as.formula(paste("p4.pray", x, " + rel3x")), data = working)
rel3x.behave <- lm(as.formula(paste("behave", x, " + rel3x")), data = working)
rel3x.click <- lm(as.formula(paste("clicked", x, " + rel3x")), data = working2)

## Mediation
## Rel Index 1, Endorsement Condition
rel1x.e.therm <- mediate(m.rel1x, rel1x.therm, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel1x.e.p1 <- mediate(m.rel1x, rel1x.p1, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel1x.e.p2 <- mediate(m.rel1x, rel1x.p2, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel1x.e.p3 <- mediate(m.rel1x, rel1x.p3, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel1x.e.p4 <- mediate(m.rel1x, rel1x.p4, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel1x.e.behave <- mediate(m.rel1x, rel1x.behave, treat = "treat", mediator = "rel1x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel1x.e.click <- mediate(m2.rel1x, rel1x.click, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## Rel Index 1, Norms Condition
rel1x.n.therm <- mediate(m.rel1x, rel1x.therm, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel1x.n.p1 <- mediate(m.rel1x, rel1x.p1, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel1x.n.p2 <- mediate(m.rel1x, rel1x.p2, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel1x.n.p3 <- mediate(m.rel1x, rel1x.p3, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel1x.n.p4 <- mediate(m.rel1x, rel1x.p4, treat = "treat", mediator = "rel1x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel1x.n.behave <- mediate(m.rel1x, rel1x.behave, treat = "treat", mediator = "rel1x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel1x.n.click <- mediate(m2.rel1x, rel1x.click, treat = "treat", mediator = "rel1x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## Rel Index 2, Endorsement Condition
rel2x.e.therm <- mediate(m.rel2x, rel2x.therm, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel2x.e.p1 <- mediate(m.rel2x, rel2x.p1, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel2x.e.p2 <- mediate(m.rel2x, rel2x.p2, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel2x.e.p3 <- mediate(m.rel2x, rel2x.p3, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel2x.e.p4 <- mediate(m.rel2x, rel2x.p4, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel2x.e.behave <- mediate(m.rel2x, rel2x.behave, treat = "treat", mediator = "rel2x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel2x.e.click <- mediate(m2.rel2x, rel2x.click, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## Rel Index 2, Norms Condition
rel2x.n.therm <- mediate(m.rel2x, rel2x.therm, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel2x.n.p1 <- mediate(m.rel2x, rel2x.p1, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel2x.n.p2 <- mediate(m.rel2x, rel2x.p2, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel2x.n.p3 <- mediate(m.rel2x, rel2x.p3, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel2x.n.p4 <- mediate(m.rel2x, rel2x.p4, treat = "treat", mediator = "rel2x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel2x.n.behave <- mediate(m.rel2x, rel2x.behave, treat = "treat", mediator = "rel2x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel2x.n.click <- mediate(m2.rel2x, rel2x.click, treat = "treat", mediator = "rel2x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## Rel Index 3, Endorsement Condition
rel3x.e.therm <- mediate(m.rel3x, rel3x.therm, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel3x.e.p1 <- mediate(m.rel3x, rel3x.p1, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel3x.e.p2 <- mediate(m.rel3x, rel3x.p2, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel3x.e.p3 <- mediate(m.rel3x, rel3x.p3, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel3x.e.p4 <- mediate(m.rel3x, rel3x.p4, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel3x.e.behave <- mediate(m.rel3x, rel3x.behave, treat = "treat", mediator = "rel3x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

rel3x.e.click <- mediate(m2.rel3x, rel3x.click, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Endorsement", robustSE = TRUE, sims = 1000)

## Rel Index 3, Norms Condition
rel3x.n.therm <- mediate(m.rel3x, rel3x.therm, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel3x.n.p1 <- mediate(m.rel3x, rel3x.p1, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel3x.n.p2 <- mediate(m.rel3x, rel3x.p2, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel3x.n.p3 <- mediate(m.rel3x, rel3x.p3, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel3x.n.p4 <- mediate(m.rel3x, rel3x.p4, treat = "treat", mediator = "rel3x", 
                      dropobs = TRUE, control.value = "Control", 
                      treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel3x.n.behave <- mediate(m.rel3x, rel3x.behave, treat = "treat", mediator = "rel3x", 
                          dropobs = TRUE, control.value = "Control", 
                          treat.value  = "Norms", robustSE = TRUE, sims = 1000)

rel3x.n.click <- mediate(m2.rel3x, rel3x.click, treat = "treat", mediator = "rel3x", 
                         dropobs = TRUE, control.value = "Control", 
                         treat.value  = "Norms", robustSE = TRUE, sims = 1000)

## Mediation Plots (rel1x)
## Rel Index 1, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel1x.e.therm, main = "Thermometer")
plot(rel1x.e.p1, main = "Policy: Hand Washing")
plot(rel1x.e.p2, main = "Policy: Distancing")
plot(rel1x.e.p3, main = "Policy: Mask Wearing")
plot(rel1x.e.p4, main = "Policy: Prayer")
plot(rel1x.e.behave, main = "Would Click")
plot(rel1x.e.click, main = "Clicked")
mtext("Mediator Variable: Regularly Attends Church \n Endorsement Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 1, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel1x.n.therm, main = "Thermometer")
plot(rel1x.n.p1, main = "Policy: Hand Washing")
plot(rel1x.n.p2, main = "Policy: Distancing")
plot(rel1x.n.p3, main = "Policy: Mask Wearing")
plot(rel1x.n.p4, main = "Policy: Prayer")
plot(rel1x.n.behave, main = "Would Click")
plot(rel1x.n.click, main = "Clicked")
mtext("Mediator Variable: Regularly Attends Church \n Norms Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 2, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel2x.e.therm, main = "Thermometer")
plot(rel2x.e.p1, main = "Policy: Hand Washing")
plot(rel2x.e.p2, main = "Policy: Distancing")
plot(rel2x.e.p3, main = "Policy: Mask Wearing")
plot(rel2x.e.p4, main = "Policy: Prayer")
plot(rel2x.e.behave, main = "Would Click")
plot(rel2x.e.click, main = "Clicked")
mtext("Mediator Variable: Spiritual Values \n Endorsement Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 2, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel2x.n.therm, main = "Thermometer")
plot(rel2x.n.p1, main = "Policy: Hand Washing")
plot(rel2x.n.p2, main = "Policy: Distancing")
plot(rel2x.n.p3, main = "Policy: Mask Wearing")
plot(rel2x.n.p4, main = "Policy: Prayer")
plot(rel2x.n.behave, main = "Would Click")
plot(rel2x.n.click, main = "Clicked")
mtext("Mediator Variable: Spiritual Values \n Norms Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 3, Endorsement Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel3x.e.therm, main = "Thermometer")
plot(rel3x.e.p1, main = "Policy: Hand Washing")
plot(rel3x.e.p2, main = "Policy: Distancing")
plot(rel3x.e.p3, main = "Policy: Mask Wearing")
plot(rel3x.e.p4, main = "Policy: Prayer")
plot(rel3x.e.behave, main = "Would Click")
plot(rel3x.e.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Endorsement Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Rel Index 3, Norms Condition
par(mfrow=c(3,3), oma=c(0,0,3,0))
plot(rel3x.n.therm, main = "Thermometer")
plot(rel3x.n.p1, main = "Policy: Hand Washing")
plot(rel3x.n.p2, main = "Policy: Distancing")
plot(rel3x.n.p3, main = "Policy: Mask Wearing")
plot(rel3x.n.p4, main = "Policy: Prayer")
plot(rel3x.n.behave, main = "Would Click")
plot(rel3x.n.click, main = "Clicked")
mtext("Mediator Variable: Enjoys Social Aspects of Church \n Norms Condition", outer=TRUE,  side = 3, cex=1, line=-0.5)

## Cleaning Up
rm(list=setdiff(ls(), c("working", "working2")))
```

## Interaction Models

We also explore the religiosity and religious motivation mechanisms by estimating conditional average treatment effects with a factorial interaction model.  We present these results below.

### Religiosity (Self Reported)

```{R int_religiosity, fig.width = 10, fig.height = 8}
## Religiosity Mechanism #######################################################
## Data Cleaning
x <- " ~ treat:religiosity + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models
int.therm <- lm(as.formula(paste("therm", x)), data = working)
int.p1 <- lm(as.formula(paste("p1.wash", x)), data = working)
int.p2 <- lm(as.formula(paste("p2.dist", x)), data = working)
int.p3 <- lm(as.formula(paste("p3.mask", x)), data = working)
int.p4 <- lm(as.formula(paste("p4.pray", x)), data = working)
int.behave <- lm(as.formula(paste("behave", x)), data = working)
int.click <- lm(as.formula(paste("clicked", x)), data = working2)

## OLS Tables
int.mod <- list("Thermometer" = int.therm, 
                "Policy: Hand Washing" = int.p1, 
                "Policy: Distancing" = int.p2, 
                "Policy: Mask Wearing" = int.p3, 
                "Policy: Prayer" = int.p4, 
                "Would Click" = int.behave, 
                "Clicked" = int.click)

cm <- c("(Intercept)" =	"Intercept",
        "treatNorms:religiosityAnti-Religious" =	"Norms : Anti-Religious",
        "treatNorms:religiosityNot at all Religious" =	"Norms : Not at all Religious",
        "treatNorms:religiositySlightly Religious" =	"Norms : Slightly Religious",
        "treatNorms:religiosityModerately Religious" =	"Norms : Moderately Religious",
        "treatNorms:religiosityVery Religious" =	"Norms : Very Religious",        
        "treatEndorsement:religiosityAnti-Religious" =	"Endorsement : Anti-Religious",
        "treatEndorsement:religiosityNot at all Religious" =	"Endorsement : Not at all Religious",
        "treatEndorsement:religiositySlightly Religious" =	"Endorsement : Slightly Religious",
        "treatEndorsement:religiosityModerately Religious" =	"Endorsement : Moderately Religious",
        "treatEndorsement:religiosityVery Religious" =	"Endorsement : Very Religious",
        "treatControl:religiosityAnti-Religious" =	"Control : Anti-Religious",
        "treatControl:religiosityNot at all Religious" =	"Control : Not at all Religious",        
        "treatControl:religiositySlightly Religious" =	"Control : Slightly Religious",
        "treatControl:religiosityModerately Religious" =	"Control : Moderately Religious",
        "treatControl:religiosityVery Religious" =	"Control : Very Religious")

modelsummary(int.mod,
             title = "Interaction Model: Treatment x Religiosity (Self-Reported),
             Omitted Category = Control x Anti-Religious",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Margins Plot
p.int.therm <- ggpredict(int.therm, terms = c("treat", "religiosity"))
p.int.p1 <- ggpredict(int.p1, terms = c("treat", "religiosity"))
p.int.p2 <- ggpredict(int.p2, terms = c("treat", "religiosity"))
p.int.p3 <- ggpredict(int.p3, terms = c("treat", "religiosity"))
p.int.p4 <- ggpredict(int.p4, terms = c("treat", "religiosity"))
p.int.behave <- ggpredict(int.behave, terms = c("treat", "religiosity"))
p.int.click <- ggpredict(int.click, terms = c("treat", "religiosity"))

p.int.therm$model <- "Thermometer"
p.int.p1$model <- "Policy: Hand Washing"
p.int.p2$model <- "Policy: Distancing"
p.int.p3$model <- "Policy: Mask Wearing"
p.int.p4$model <- "Policy: Prayer"
p.int.behave$model <- "Would Click"
p.int.click$model <- "Clicked"

p.int.all <- rbind(p.int.therm, p.int.p1, p.int.p2, p.int.p3, p.int.p4, p.int.behave, p.int.click)

p.int.all$model <- factor(p.int.all$model, 
                          levels = c("Thermometer", 
                                     "Policy: Hand Washing", 
                                     "Policy: Distancing", 
                                     "Policy: Mask Wearing", 
                                     "Policy: Prayer", 
                                     "Would Click", 
                                     "Clicked"), 
                          ordered = TRUE)

p.int.all$x <- factor(p.int.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

ggplot(p.int.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Religiosity (Self-Reported)",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), c("working")))
```

### Religious Motivations

```{R int_rel_motiv, fig.width = 10, fig.height = 8}

## Religious Motivations #######################################################
## Create a factor variable for the 6 religious motivations
# as.data.frame(colnames(working))

working2 <- working
working2$m6 <- factor(apply(working2[,51:56], 1, function(x) names(x)[x == 1]), 
                      levels = names(working2[,51:56]))

working3 <- working2 %>%
  filter(behave > 4)

working2$treat <- factor(working2$treat, levels = c("Norms", "Endorsement", "Control"))
working3$treat <- factor(working3$treat, levels = c("Norms", "Endorsement", "Control"))

## Note: Tables XXX in SI-F present the same results but changed the...
## ...omitted (reference) category
## To change the reference category for the interaction term:
# working2$treat <- relevel(working2$treat, ref = "Control") # Change Category Here
# working3$treat <- relevel(working3$treat, ref = "Control") # Change Category Here
# working2$m6 <- relevel(working2$m6, ref = "c1") # Change Category Here
# working3$m6 <- relevel(working3$m6, ref = "c1") # Change Category Here

x <- "~ m6:treat + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models
int.therm <- lm(as.formula(paste("therm", x)), data = working2)
int.p1 <- lm(as.formula(paste("p1.wash", x)), data = working2)
int.p2 <- lm(as.formula(paste("p2.dist", x)), data = working2)
int.p3 <- lm(as.formula(paste("p3.mask", x)), data = working2)
int.p4 <- lm(as.formula(paste("p4.pray", x)), data = working2)
int.behave <- lm(as.formula(paste("behave", x)), data = working2)
int.click <- lm(as.formula(paste("clicked", x)), data = working3)

## OLS Tables
m6.mod <- list("Thermometer" = int.therm, 
               "Policy: Hand Washing" = int.p1, 
               "Policy: Distancing" = int.p2, 
               "Policy: Mask Wearing" = int.p3, 
               "Policy: Prayer" = int.p4, 
               "Would Click" = int.behave, 
               "Clicked" = int.click)

cm <- c("(Intercept)" = "Intercept",
        "m6c1:treatNorms" = "Norms : Insurance",
        "m6c2:treatNorms" = "Norms : Morality",
        "m6c3:treatNorms" = "Norms : Social Capital",
        "m6c4:treatNorms" = "Norms : Heuristics",
        "m6c5:treatNorms" = "Norms : Tradition",
        "m6c6:treatNorms" = "Norms : Social Identity",
        "m6c1:treatEndorsement" = "Endorsement : Insurance",
        "m6c2:treatEndorsement" = "Endorsement : Morality",
        "m6c3:treatEndorsement" = "Endorsement : Social Capital",
        "m6c4:treatEndorsement" = "Endorsement : Heuristics",
        "m6c5:treatEndorsement" = "Endorsement : Tradition",
        "m6c6:treatEndorsement" = "Endorsement : Social Identity",
        "m6c1:treatControl" = "Control : Insurance",
        "m6c2:treatControl" = "Control : Morality",
        "m6c3:treatControl" = "Control : Social Capital",
        "m6c4:treatControl" = "Control : Heuristics",
        "m6c5:treatControl" = "Control : Tradition",
        "m6c6:treatControl" = "Control : Social Identity")

modelsummary(m6.mod,
             title = "Interaction Model: Treatment x Religious Motivations, 
             Omitted Category = Control x Social Identity",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Margins Plot
m.therm <- ggpredict(int.therm, terms = c("treat", "m6"))
m.therm$model <- "Thermometer"

m.p1 <- ggpredict(int.p1, terms = c("treat", "m6"))
m.p1$model <- "Policy: Hand Washing"

m.p2 <- ggpredict(int.p2, terms = c("treat", "m6"))
m.p2$model <- "Policy: Distancing"

m.p3 <- ggpredict(int.p3, terms = c("treat", "m6"))
m.p3$model <- "Policy: Mask Wearing"

m.p4 <- ggpredict(int.p4, terms = c("treat", "m6"))
m.p4$model <- "Policy: Prayer"

m.behave <- ggpredict(int.behave, terms = c("treat", "m6"))
m.behave$model <- "Would Click"

m.click <- ggpredict(int.click, terms = c("treat", "m6"))
m.click$model <- "Clicked"

m.combined <- rbind(m.therm, m.p1, m.p2, m.p3, m.p4, m.behave, m.click)

m.combined$model <- factor(m.combined$model, 
                           levels = c("Thermometer", 
                                      "Policy: Hand Washing", 
                                      "Policy: Distancing", 
                                      "Policy: Mask Wearing", 
                                      "Policy: Prayer", 
                                      "Would Click", 
                                      "Clicked"), 
                           ordered = TRUE)

m.combined$group <- dplyr::recode(m.combined$group,
                                  "c1" = "Insurance",
                                  "c2" = "Morality",
                                  "c3" = "Social Capital",
                                  "c4" = "Heuristics",
                                  "c5" = "Tradition",
                                  "c6" = "Social Identity")

ggplot(m.combined, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.5, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Religious Motivations",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), c("working", "working2")))
```

### Extrinsic / Intrinsic Motivations

```{R int_ei_motiv, fig.width = 10, fig.height = 8}
## Extrinsic / Intrinsic Motivation ############################################
working2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat:ei1.f + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat:ei2.f + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv3 <- " ~ treat:ei3.f + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models & Tables - EI 1 - Religion Offers Comfort
i1.therm <- lm(as.formula(paste("therm", iv1)), data = working)
i1.p1 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
i1.p2 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
i1.p3 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
i1.p4 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
i1.behave <- lm(as.formula(paste("behave", iv1)), data = working)
i1.click <- lm(as.formula(paste("clicked", iv1)), data = working2)

i1 <- list("Thermometer" = i1.therm, 
           "Policy: Hand Washing" = i1.p1, 
           "Policy: Distancing" = i1.p2, 
           "Policy: Mask Wearing" = i1.p3, 
           "Policy: Prayer" = i1.p4, 
           "Would Click" = i1.behave, 
           "Clicked" = i1.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:ei1.fStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:ei1.fDisagree" = "Norms : Disagree",
        "treatNorms:ei1.fNeither Agree nor Disagree" = "Norms : Neither Agree nor Disagree",
        "treatNorms:ei1.fAgree" = "Norms : Agree",
        "treatNorms:ei1.fStrongly Agree" = "Norms : Stongly Agree",
        "treatEndorsement:ei1.fStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:ei1.fDisagree" = "Endorsement : Disagree",
        "treatEndorsement:ei1.fNeither Agree nor Disagree" = "Endorsement : Neither Agree nor Disagree",
        "treatEndorsement:ei1.fAgree" = "Endorsement : Agree",
        "treatEndorsement:ei1.fStrongly Agree" = "Endorsement : Stongly Agree",
      	"treatControl:ei1.fStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:ei1.fDisagree" = "Control : Disagree",        
        "treatControl:ei1.fNeither Agree nor Disagree" = "Control : Neither Agree nor Disagree",
        "treatControl:ei1.fAgree" = "Control : Agree",
        "treatControl:ei1.fStrongly Agree" = "Control : Stongly Agree")

modelsummary(i1,
             title = "Interaction Model: Treatment x EI Motivations (Religion Offers Comfort), Omitted Cateogry = Norms x Strongly Agree",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - EI 2 - Religious Approach to Life
i2.therm <- lm(as.formula(paste("therm", iv2)), data = working)
i2.p1 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
i2.p2 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
i2.p3 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
i2.p4 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
i2.behave <- lm(as.formula(paste("behave", iv2)), data = working)
i2.click <- lm(as.formula(paste("clicked", iv2)), data = working2)

i2 <- list("Thermometer" = i2.therm, 
           "Policy: Hand Washing" = i2.p1, 
           "Policy: Distancing" = i2.p2, 
           "Policy: Mask Wearing" = i2.p3, 
           "Policy: Prayer" = i2.p4, 
           "Would Click" = i2.behave, 
           "Clicked" = i2.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:ei2.fStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:ei2.fDisagree" = "Norms : Disagree",
        "treatNorms:ei2.fNeither Agree nor Disagree" = "Norms : Neither Agree nor Disagree",
        "treatNorms:ei2.fAgree" = "Norms : Agree",
        "treatNorms:ei2.fStrongly Agree" = "Norms : Stongly Agree",
        "treatEndorsement:ei2.fStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:ei2.fDisagree" = "Endorsement : Disagree",
        "treatEndorsement:ei2.fNeither Agree nor Disagree" = "Endorsement : Neither Agree nor Disagree",
        "treatEndorsement:ei2.fAgree" = "Endorsement : Agree",
        "treatEndorsement:ei2.fStrongly Agree" = "Endorsement : Stongly Agree",
    	  "treatControl:ei2.fStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:ei2.fDisagree" = "Control : Disagree",        
        "treatControl:ei2.fNeither Agree nor Disagree" = "Control : Neither Agree nor Disagree",
        "treatControl:ei2.fAgree" = "Control : Agree",
        "treatControl:ei2.fStrongly Agree" = "Control : Stongly Agree")

modelsummary(i2,
             title = "Interaction Model: Treatment x EI Motivations (Takes Religious Approach to Life), Omitted Category = Norms x Strongly Agree",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - EI 3 - Enjoys Social Aspects of Church
i3.therm <- lm(as.formula(paste("therm", iv3)), data = working)
i3.p1 <- lm(as.formula(paste("p1.wash", iv3)), data = working)
i3.p2 <- lm(as.formula(paste("p2.dist", iv3)), data = working)
i3.p3 <- lm(as.formula(paste("p3.mask", iv3)), data = working)
i3.p4 <- lm(as.formula(paste("p4.pray", iv3)), data = working)
i3.behave <- lm(as.formula(paste("behave", iv3)), data = working)
i3.click <- lm(as.formula(paste("clicked", iv3)), data = working2)

i3 <- list("Thermometer" = i3.therm, 
           "Policy: Hand Washing" = i3.p1, 
           "Policy: Distancing" = i3.p2, 
           "Policy: Mask Wearing" = i3.p3, 
           "Policy: Prayer" = i3.p4, 
           "Would Click" = i3.behave, 
           "Clicked" = i3.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:ei3.fStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:ei3.fDisagree" = "Norms : Disagree",
        "treatNorms:ei3.fNeither Agree nor Disagree" = "Norms : Neither Agree nor Disagree",
        "treatNorms:ei3.fAgree" = "Norms : Agree",
        "treatNorms:ei3.fStrongly Agree" = "Norms : Stongly Agree",
        "treatEndorsement:ei3.fStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:ei3.fDisagree" = "Endorsement : Disagree",
        "treatEndorsement:ei3.fNeither Agree nor Disagree" = "Endorsement : Neither Agree nor Disagree",
        "treatEndorsement:ei3.fAgree" = "Endorsement : Agree",
        "treatEndorsement:ei3.fStrongly Agree" = "Endorsement : Stongly Agree",
    	  "treatControl:ei3.fStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:ei3.fDisagree" = "Control : Disagree",        
        "treatControl:ei3.fNeither Agree nor Disagree" = "Control : Neither Agree nor Disagree",
        "treatControl:ei3.fAgree" = "Control : Agree",
        "treatControl:ei3.fStrongly Agree" = "Control : Stongly Agree")

modelsummary(i3,
             title = "Interaction Model: Treatment x EI Motivations (Enjoys Social Aspects of Church), Omitted Category = Norms x Strongly Agree",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Predictive Margins for EI Interaction Models
p.i1.therm <- ggpredict(i1.therm, terms = c("treat", "ei1.f"))
p.i1.p1 <- ggpredict(i1.p1, terms = c("treat", "ei1.f"))
p.i1.p2 <- ggpredict(i1.p2, terms = c("treat", "ei1.f"))
p.i1.p3 <- ggpredict(i1.p3, terms = c("treat", "ei1.f"))
p.i1.p4 <- ggpredict(i1.p4, terms = c("treat", "ei1.f"))
p.i1.behave <- ggpredict(i1.behave, terms = c("treat", "ei1.f"))
p.i1.click <- ggpredict(i1.click, terms = c("treat", "ei1.f"))

p.i1.therm$model <- "Thermometer"
p.i1.p1$model <- "Policy: Hand Washing"
p.i1.p2$model <- "Policy: Distancing"
p.i1.p3$model <- "Policy: Mask Wearing"
p.i1.p4$model <- "Policy: Prayer"
p.i1.behave$model <- "Would Click"
p.i1.click$model <- "Clicked"

p.i1.all <- rbind(p.i1.therm, p.i1.p1, p.i1.p2, p.i1.p3, p.i1.p4, p.i1.behave, p.i1.click)

p.i1.all$x <- factor(p.i1.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

p.i1.all$model <- factor(p.i1.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i2.therm <- ggpredict(i2.therm, terms = c("treat", "ei2.f"))
p.i2.p1 <- ggpredict(i2.p1, terms = c("treat", "ei2.f"))
p.i2.p2 <- ggpredict(i2.p2, terms = c("treat", "ei2.f"))
p.i2.p3 <- ggpredict(i2.p3, terms = c("treat", "ei2.f"))
p.i2.p4 <- ggpredict(i2.p4, terms = c("treat", "ei2.f"))
p.i2.behave <- ggpredict(i2.behave, terms = c("treat", "ei2.f"))
p.i2.click <- ggpredict(i2.click, terms = c("treat", "ei2.f"))

p.i2.therm$model <- "Thermometer"
p.i2.p1$model <- "Policy: Hand Washing"
p.i2.p2$model <- "Policy: Distancing"
p.i2.p3$model <- "Policy: Mask Wearing"
p.i2.p4$model <- "Policy: Prayer"
p.i2.behave$model <- "Would Click"
p.i2.click$model <- "Clicked"

p.i2.all <- rbind(p.i2.therm, p.i2.p1, p.i2.p2, p.i2.p3, p.i2.p4, p.i2.behave, p.i2.click)

p.i2.all$x <- factor(p.i2.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

p.i2.all$model <- factor(p.i2.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i3.therm <- ggpredict(i3.therm, terms = c("treat", "ei3.f"))
p.i3.p1 <- ggpredict(i3.p1, terms = c("treat", "ei3.f"))
p.i3.p2 <- ggpredict(i3.p2, terms = c("treat", "ei3.f"))
p.i3.p3 <- ggpredict(i3.p3, terms = c("treat", "ei3.f"))
p.i3.p4 <- ggpredict(i3.p4, terms = c("treat", "ei3.f"))
p.i3.behave <- ggpredict(i3.behave, terms = c("treat", "ei3.f"))
p.i3.click <- ggpredict(i3.click, terms = c("treat", "ei3.f"))

p.i3.therm$model <- "Thermometer"
p.i3.p1$model <- "Policy: Hand Washing"
p.i3.p2$model <- "Policy: Distancing"
p.i3.p3$model <- "Policy: Mask Wearing"
p.i3.p4$model <- "Policy: Prayer"
p.i3.behave$model <- "Would Click"
p.i3.click$model <- "Clicked"

p.i3.all <- rbind(p.i3.therm, p.i3.p1, p.i3.p2, p.i3.p3, p.i3.p4, p.i3.behave, p.i3.click)

p.i3.all$x <- factor(p.i3.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

p.i3.all$model <- factor(p.i3.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

## Margins Plots
ggplot(p.i1.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Extrinsic / Intrinsic Motivations - Religion Offers Comfort",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i2.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Extrinsic / Intrinsic Motivations - Takes Religious Approach to Life",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i3.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Extrinsic / Intrinsic Motivations - Enjoys Social Aspects of Church",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), c("working", "working2")))
```

### Religiosity Index

```{R int_relindex_motiv, fig.width = 10, fig.height = 8}
## Religiosity Index ###########################################################
iv1 <- " ~ treat:rel1x + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat:rel2x + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv3 <- " ~ treat:rel3x + male + age + education + party.f + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"

## Interaction Models & Tables - Rel 1 - Regularly Attends Church
i1.therm <- lm(as.formula(paste("therm", iv1)), data = working)
i1.p1 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
i1.p2 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
i1.p3 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
i1.p4 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
i1.behave <- lm(as.formula(paste("behave", iv1)), data = working)
i1.click <- lm(as.formula(paste("clicked", iv1)), data = working2)

i1 <- list("Thermometer" = i1.therm, 
           "Policy: Hand Washing" = i1.p1, 
           "Policy: Distancing" = i1.p2, 
           "Policy: Mask Wearing" = i1.p3, 
           "Policy: Prayer" = i1.p4, 
           "Would Click" = i1.behave, 
           "Clicked" = i1.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:rel1xStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:rel1xDisagree" = "Norms : Disagree",
        "treatNorms:rel1xSomewhat Disagree" = "Norms : Somewhat Disagree",
        "treatNorms:rel1xSomewhat Agree" = "Norms : Somewhat Agree",
        "treatNorms:rel1xAgree" = "Norms : Agree",
        "treatNorms:rel1xStrongly Agree" = "Norms : Strongly Agree",
        "treatEndorsement:rel1xStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:rel1xDisagree" = "Endorsement : Disagree",
        "treatEndorsement:rel1xSomewhat Disagree" = "Endorsement : Somewhat Disagree",
        "treatEndorsement:rel1xSomewhat Agree" = "Endorsement : Somewhat Agree",
        "treatEndorsement:rel1xAgree" = "Endorsement : Agree",
        "treatEndorsement:rel1xStrongly Agree" = "Endorsement : Strongly Agree",
        "treatControl:rel1xStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:rel1xDisagree" = "Control : Disagree",
        "treatControl:rel1xSomewhat Disagree" = "Control : Somewhat Disagree",
        "treatControl:rel1xSomewhat Agree" = "Control : Somewhat Agree",
        "treatControl:rel1xAgree" = "Control : Agree",
        "treatControl:rel1xStrongly Agree" = "Control : Strongly Agree")

modelsummary(i1,
             title = "Interaction Model: Treatment x Regularly Attends Church, Omitted Category = Norms x Strongly Agree",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - Rel 2 - Spiritual Values
i2.therm <- lm(as.formula(paste("therm", iv2)), data = working)
i2.p1 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
i2.p2 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
i2.p3 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
i2.p4 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
i2.behave <- lm(as.formula(paste("behave", iv2)), data = working)
i2.click <- lm(as.formula(paste("clicked", iv2)), data = working2)

i2 <- list("Thermometer" = i2.therm, 
           "Policy: Hand Washing" = i2.p1, 
           "Policy: Distancing" = i2.p2, 
           "Policy: Mask Wearing" = i2.p3, 
           "Policy: Prayer" = i2.p4, 
           "Would Click" = i2.behave, 
           "Clicked" = i2.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:rel2xStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:rel2xDisagree" = "Norms : Disagree",
        "treatNorms:rel2xSomewhat Disagree" = "Norms : Somewhat Disagree",
        "treatNorms:rel2xSomewhat Agree" = "Norms : Somewhat Agree",
        "treatNorms:rel2xAgree" = "Norms : Agree",
        "treatNorms:rel2xStrongly Agree" = "Norms : Strongly Agree",
        "treatEndorsement:rel2xStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:rel2xDisagree" = "Endorsement : Disagree",
        "treatEndorsement:rel2xSomewhat Disagree" = "Endorsement : Somewhat Disagree",
        "treatEndorsement:rel2xSomewhat Agree" = "Endorsement : Somewhat Agree",
        "treatEndorsement:rel2xAgree" = "Endorsement : Agree",
        "treatEndorsement:rel2xStrongly Agree" = "Endorsement : Strongly Agree",
        "treatControl:rel2xStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:rel2xDisagree" = "Control : Disagree",
        "treatControl:rel2xSomewhat Disagree" = "Control : Somewhat Disagree",
        "treatControl:rel2xSomewhat Agree" = "Control : Somewhat Agree",
        "treatControl:rel2xAgree" = "Control : Agree",
        "treatControl:rel2xStrongly Agree" = "Control : Strongly Agree")

modelsummary(i2,
             title = "Interaction Model: Treatment x Spiritual Values, Omitted Category = Norms x Strongly Agree",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Interaction Models & Tables - Rel 3 - Americans Should be More Religious
i3.therm <- lm(as.formula(paste("therm", iv3)), data = working)
i3.p1 <- lm(as.formula(paste("p1.wash", iv3)), data = working)
i3.p2 <- lm(as.formula(paste("p2.dist", iv3)), data = working)
i3.p3 <- lm(as.formula(paste("p3.mask", iv3)), data = working)
i3.p4 <- lm(as.formula(paste("p4.pray", iv3)), data = working)
i3.behave <- lm(as.formula(paste("behave", iv3)), data = working)
i3.click <- lm(as.formula(paste("clicked", iv3)), data = working2)

i3 <- list("Thermometer" = i3.therm, 
           "Policy: Hand Washing" = i3.p1, 
           "Policy: Distancing" = i3.p2, 
           "Policy: Mask Wearing" = i3.p3, 
           "Policy: Prayer" = i3.p4, 
           "Would Click" = i3.behave, 
           "Clicked" = i3.click)

cm <- c("(Intercept)" = "Intercept",
        "treatNorms:rel3xStrongly Disagree" = "Norms : Strongly Disagree",
        "treatNorms:rel3xDisagree" = "Norms : Disagree",
        "treatNorms:rel3xSomewhat Disagree" = "Norms : Somewhat Disagree",
        "treatNorms:rel3xSomewhat Agree" = "Norms : Somewhat Agree",
        "treatNorms:rel3xAgree" = "Norms : Agree",
        "treatNorms:rel3xStrongly Agree" = "Norms : Strongly Agree",
        "treatEndorsement:rel3xStrongly Disagree" = "Endorsement : Strongly Disagree",
        "treatEndorsement:rel3xDisagree" = "Endorsement : Disagree",
        "treatEndorsement:rel3xSomewhat Disagree" = "Endorsement : Somewhat Disagree",
        "treatEndorsement:rel3xSomewhat Agree" = "Endorsement : Somewhat Agree",
        "treatEndorsement:rel3xAgree" = "Endorsement : Agree",
        "treatEndorsement:rel3xStrongly Agree" = "Endorsement : Strongly Agree",
        "treatControl:rel3xStrongly Disagree" = "Control : Strongly Disagree",
        "treatControl:rel3xDisagree" = "Control : Disagree",
        "treatControl:rel3xSomewhat Disagree" = "Control : Somewhat Disagree",
        "treatControl:rel3xSomewhat Agree" = "Control : Somewhat Agree",
        "treatControl:rel3xAgree" = "Control : Agree",
        "treatControl:rel3xStrongly Agree" = "Control : Strongly Agree")
        
modelsummary(i3,
             title = "Interaction Model: Treatment × Americans Should be More Religious, Omitted Category = Norms x Strongly Agree",
             stars = TRUE,
             coef_map = cm,
             gof_omit = 'IC|Log|Adj')

## Predictive Margins for Religiosity Index Interaction Models
p.i1.therm <- ggpredict(i1.therm, terms = c("treat", "rel1x"))
p.i1.p1 <- ggpredict(i1.p1, terms = c("treat", "rel1x"))
p.i1.p2 <- ggpredict(i1.p2, terms = c("treat", "rel1x"))
p.i1.p3 <- ggpredict(i1.p3, terms = c("treat", "rel1x"))
p.i1.p4 <- ggpredict(i1.p4, terms = c("treat", "rel1x"))
p.i1.behave <- ggpredict(i1.behave, terms = c("treat", "rel1x"))
p.i1.click <- ggpredict(i1.click, terms = c("treat", "rel1x"))

p.i1.therm$model <- "Thermometer"
p.i1.p1$model <- "Policy: Hand Washing"
p.i1.p2$model <- "Policy: Distancing"
p.i1.p3$model <- "Policy: Mask Wearing"
p.i1.p4$model <- "Policy: Prayer"
p.i1.behave$model <- "Would Click"
p.i1.click$model <- "Clicked"

p.i1.all <- rbind(p.i1.therm, p.i1.p1, p.i1.p2, p.i1.p3, p.i1.p4, p.i1.behave, p.i1.click)

p.i1.all$x <- factor(p.i1.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

p.i1.all$model <- factor(p.i1.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i2.therm <- ggpredict(i2.therm, terms = c("treat", "rel2x"))
p.i2.p1 <- ggpredict(i2.p1, terms = c("treat", "rel2x"))
p.i2.p2 <- ggpredict(i2.p2, terms = c("treat", "rel2x"))
p.i2.p3 <- ggpredict(i2.p3, terms = c("treat", "rel2x"))
p.i2.p4 <- ggpredict(i2.p4, terms = c("treat", "rel2x"))
p.i2.behave <- ggpredict(i2.behave, terms = c("treat", "rel2x"))
p.i2.click <- ggpredict(i2.click, terms = c("treat", "rel2x"))

p.i2.therm$model <- "Thermometer"
p.i2.p1$model <- "Policy: Hand Washing"
p.i2.p2$model <- "Policy: Distancing"
p.i2.p3$model <- "Policy: Mask Wearing"
p.i2.p4$model <- "Policy: Prayer"
p.i2.behave$model <- "Would Click"
p.i2.click$model <- "Clicked"

p.i2.all <- rbind(p.i2.therm, p.i2.p1, p.i2.p2, p.i2.p3, p.i2.p4, p.i2.behave, p.i2.click)

p.i2.all$x <- factor(p.i2.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

p.i2.all$model <- factor(p.i2.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)

p.i3.therm <- ggpredict(i3.therm, terms = c("treat", "rel3x"))
p.i3.p1 <- ggpredict(i3.p1, terms = c("treat", "rel3x"))
p.i3.p2 <- ggpredict(i3.p2, terms = c("treat", "rel3x"))
p.i3.p3 <- ggpredict(i3.p3, terms = c("treat", "rel3x"))
p.i3.p4 <- ggpredict(i3.p4, terms = c("treat", "rel3x"))
p.i3.behave <- ggpredict(i3.behave, terms = c("treat", "rel3x"))
p.i3.click <- ggpredict(i3.click, terms = c("treat", "rel3x"))

p.i3.therm$model <- "Thermometer"
p.i3.p1$model <- "Policy: Hand Washing"
p.i3.p2$model <- "Policy: Distancing"
p.i3.p3$model <- "Policy: Mask Wearing"
p.i3.p4$model <- "Policy: Prayer"
p.i3.behave$model <- "Would Click"
p.i3.click$model <- "Clicked"

p.i3.all <- rbind(p.i3.therm, p.i3.p1, p.i3.p2, p.i3.p3, p.i3.p4, p.i3.behave, p.i3.click)

p.i3.all$x <- factor(p.i3.all$x, 
                levels = c("Norms", 
                           "Endorsement", 
                           "Control") ,
                ordered = TRUE)

p.i3.all$model <- factor(p.i3.all$model, 
                         levels = c("Thermometer", 
                                    "Policy: Hand Washing", 
                                    "Policy: Distancing", 
                                    "Policy: Mask Wearing", 
                                    "Policy: Prayer", 
                                    "Would Click", 
                                    "Clicked"), 
                         ordered = TRUE)


## Plots
ggplot(p.i1.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Regularly Attends Church",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i2.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Spiritual Values",
       x  = "Condition", y = "Predicted Value")

ggplot(p.i3.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model)) +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  theme(strip.text.y = element_text(size = 7)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Americans Should be More Religious",
       x  = "Condition", y = "Predicted Value")

## Cleaning Up
rm(list=setdiff(ls(), "working"))
```

## Exploratory Analysis of Partisanship

Although we did not pre-register any specifications that explicitly consider partisan identification in our pre-analysis plan, we perform an exploratory analysis in the revised draft to examine the effects of partisanship.  Specifically, we estimate the conditional average treatment effect (CATE) by party ID using the following interaction model: 

$${\textrm{Outcome}} = \alpha + \beta_1\textrm{Treat} + \beta_2\textrm{Party ID} + \beta_3(\textrm{Treat}\times\textrm{Party ID}) + \epsilon_i$$
Where *Treat* is assignment to treatment and *Party ID* is the respondent's party ID.  For ease of interpretation, we discuss a specification that uses a collapsed, three-category party variable in the manuscript, though we present results from specifications that use both the three-category and the seven-category coding in the Supplementary Information document.  

### Three Category Party ID

```{R party_3cat_1, fig.width = 10, fig.height = 8}
# Exploratory Analysis of Partisanship #########################################
## Conditional Average Treatment Effects #######################################
# levels(working$party.f)   # Change the reference category to independent
working$party.f2 <- factor(working$party.f, levels = c("Strong Democrat", "Democrat", 
                                                       "Lean Democrat", "Independent",
                                                       "Lean Republican", 
                                                       "Republican", "Strong Republican"))
working$party.32 <- factor(working$party.3, levels = c("Democrat",
                                                       "Independent", 
                                                       "Republican"))

w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat * party.32 + gender + age + education + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat * party.f2 + gender + age + education + news.source + news.follow + covid.source + covid.follow + covid.threat + covid.tested + covid.friends"


### Three Category Party Analysis ##############################################
## Simplified Model, No Controls
# a1 <- lm(therm ~ treat * party.32, data = working)
# a2 <- lm(p1.wash ~ treat * party.32, data = working)
# a3 <- lm(p2.dist ~ treat * party.32, data = working)
# a4 <- lm(p3.mask ~ treat * party.32, data = working)
# a5 <- lm(p4.pray ~ treat * party.32, data = working)
# a6 <- lm(behave ~ treat * party.32, data = working)
# a7 <- lm(clicked ~ treat * party.32, data = w2)

## Model With All Controls 
a1 <- lm(as.formula(paste("therm", iv1)), data = working)
a2 <- lm(as.formula(paste("p1.wash", iv1)), data = working)
a3 <- lm(as.formula(paste("p2.dist", iv1)), data = working)
a4 <- lm(as.formula(paste("p3.mask", iv1)), data = working)
a5 <- lm(as.formula(paste("p4.pray", iv1)), data = working)
a6 <- lm(as.formula(paste("behave", iv1)), data = working)
a7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

### 3 Cat Margins & Margins Plots
int.therm <- ggpredict(a1, terms = c("treat", "party.32"))
int.p1 <- ggpredict(a2, terms = c("treat", "party.32"))
int.p2 <- ggpredict(a3, terms = c("treat", "party.32"))
int.p3 <- ggpredict(a4, terms = c("treat", "party.32"))
int.p4 <- ggpredict(a5, terms = c("treat", "party.32"))
int.behave <- ggpredict(a6, terms = c("treat", "party.32"))
int.click <- ggpredict(a7, terms = c("treat", "party.32"))

int.therm$model <- "Thermometer"
int.p1$model <- "Policy: Hand Washing"
int.p2$model <- "Policy: Distancing"
int.p3$model <- "Policy: Mask Wearing"
int.p4$model <- "Policy: Prayer"
int.behave$model <- "Would Click"
int.click$model <- "Clicked"

int.all <- rbind(int.therm, int.p1, int.p2, int.p3, int.p4, int.behave, int.click)

int.all$model <- factor(int.all$model, 
                        levels = c("Thermometer", 
                                   "Policy: Hand Washing", 
                                   "Policy: Distancing", 
                                   "Policy: Mask Wearing", 
                                   "Policy: Prayer", 
                                   "Would Click", 
                                   "Clicked"), 
                        ordered = TRUE)

ggplot(int.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model),
             scales = "free_y") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Party (Three Category)",
       x  = "Condition", y = "Predicted Value",
       caption = "Interaction model with controls.")

rm(int.all, int.therm, int.p1, int.p2, int.p3, int.p4, int.behave, int.click)
```

The following contrast (pairwise $p$ value) plots depict the estimated effect size by party ID, along with the $p$ value associated with that comparison.

```{R party_3cat_2, fig.width = 10}
## 3 Cat Contrasts & PWPPs 
## Use this for models without controls:
# z.a1 <- emmeans(a1, ~ treat*party.32)
# z.a2 <- emmeans(a2, ~ treat*party.32)
# z.a3 <- emmeans(a3, ~ treat*party.32)
# z.a4 <- emmeans(a4, ~ treat*party.32)
# z.a5 <- emmeans(a5, ~ treat*party.32)
# z.a6 <- emmeans(a6, ~ treat*party.32)
# z.a7 <- emmeans(a7, ~ treat*party.32)

## Use this for models with controls:
z.a1 <- emmeans(a1, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a2 <- emmeans(a2, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a3 <- emmeans(a3, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a4 <- emmeans(a4, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a5 <- emmeans(a5, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a6 <- emmeans(a6, ~ treat*party.32, non.nuis = c("treat", "party.32"))
z.a7 <- emmeans(a7, ~ treat*party.32, non.nuis = c("treat", "party.32"))

pwpp(z.a1, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Thermometer") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a2, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Hand Washing") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a3, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Distancing") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a4, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Mask Wearing") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a5, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Policy: Prayer") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a6, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Would Click") +
  ggplot2::facet_grid(~ party.32)

pwpp(z.a7, by = "party.32", 
     xlab = "Tukey-Adjusted P Value",
     ylab = "Estimated Marginal Mean") +
  ggplot2::theme_bw() +
  ggplot2::ggtitle("Pairwise Comparisons", 
                   subtitle = "Clicked") +
  ggplot2::facet_grid(~ party.32)
```

### Seven Cateogry Party ID

```{R party_7cat, fig.width = 10, fig.height = 10}
### Seven Category Party Analysis ##############################################
## Simplified Model, No Controls
# a1 <- lm(therm ~ treat * party.f2, data = working)
# a2 <- lm(p1.wash ~ treat * party.f2, data = working)
# a3 <- lm(p2.dist ~ treat * party.f2, data = working)
# a4 <- lm(p3.mask ~ treat * party.f2, data = working)
# a5 <- lm(p4.pray ~ treat * party.f2, data = working)
# a6 <- lm(behave ~ treat * party.f2, data = working)
# a7 <- lm(clicked ~ treat * party.f2, data = w2)

## Model With All Controls
a1 <- lm(as.formula(paste("therm", iv2)), data = working)
a2 <- lm(as.formula(paste("p1.wash", iv2)), data = working)
a3 <- lm(as.formula(paste("p2.dist", iv2)), data = working)
a4 <- lm(as.formula(paste("p3.mask", iv2)), data = working)
a5 <- lm(as.formula(paste("p4.pray", iv2)), data = working)
a6 <- lm(as.formula(paste("behave", iv2)), data = working)
a7 <- lm(as.formula(paste("clicked", iv2)), data = w2)

## 7 Cat Margins & Margins Plots
int.therm <- ggpredict(a1, terms = c("treat", "party.f2"))
int.p1 <- ggpredict(a2, terms = c("treat", "party.f2"))
int.p2 <- ggpredict(a3, terms = c("treat", "party.f2"))
int.p3 <- ggpredict(a4, terms = c("treat", "party.f2"))
int.p4 <- ggpredict(a5, terms = c("treat", "party.f2"))
int.behave <- ggpredict(a6, terms = c("treat", "party.f2"))
int.click <- ggpredict(a7, terms = c("treat", "party.f2"))

int.therm$model <- "Thermometer"
int.p1$model <- "Policy: Hand Washing"
int.p2$model <- "Policy: Distancing"
int.p3$model <- "Policy: Mask Wearing"
int.p4$model <- "Policy: Prayer"
int.behave$model <- "Would Click"
int.click$model <- "Clicked"

int.all <- rbind(int.therm, int.p1, int.p2, int.p3, int.p4, int.behave, int.click)

int.all$model <- factor(int.all$model, 
                        levels = c("Thermometer", 
                                   "Policy: Hand Washing", 
                                   "Policy: Distancing", 
                                   "Policy: Mask Wearing", 
                                   "Policy: Prayer", 
                                   "Would Click", 
                                   "Clicked"), 
                        ordered = TRUE)

ggplot(int.all, aes(x, predicted)) +
  geom_point() +
  geom_errorbar(aes(ymin=conf.low, ymax=conf.high), width=0.1) +
  geom_text(aes(label=round(predicted, digits = 2)), 
            hjust = 1.55, vjust = 0.5, size = 1.75) +
  facet_grid(rows = vars(group), 
             cols = vars(model),
             scales = "free_y") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank()) +
  theme(axis.text.x = element_text(angle = 45, hjust = 1)) +
  labs(title = "Marginal Effects: Interaction Model",
       subtitle = "Treatment :: Party (Seven Category)",
       x  = "Condition", y = "Predicted Value",
       caption = "Interaction model with controls.") +
  theme(strip.text.y = element_text(size = 7))

## Cleaning Up
rm(list=setdiff(ls(), "working"))
```

## Republican Subgroup Analysis

To further explore the influence of partisanship, we restrict our sample to only those respondents that self-identify as Republican (lean Republican, Republican, and strong Republican) and re-run our ATE analysis. 

```{R subgroup, fig.width = 10, fig.height = 8}
## Subgroup Analysis ###########################################################
w <- working %>% 
  filter(party.3 == "Republican")

w2 <- working %>% 
  filter(behave > 4)

iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"


### Baseline OLS Models ########################################################
## The Models
a1 <- lm(therm ~ treat, data = w)
a2 <- lm(p1.wash ~ treat, data = w)
a3 <- lm(p2.dist ~ treat, data = w)
a4 <- lm(p3.mask ~ treat, data = w)
a5 <- lm(p4.pray ~ treat, data = w)
a6 <- lm(behave ~ treat, data = working)
a7 <- lm(clicked ~ treat, data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(a1, conf.int = TRUE)
x2 <- tidy(a2, conf.int = TRUE)
x3 <- tidy(a3, conf.int = TRUE)
x4 <- tidy(a4, conf.int = TRUE)
x5 <- tidy(a5, conf.int = TRUE)
x6 <- tidy(a6, conf.int = TRUE)
x7 <- tidy(a7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(a1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(a2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(a3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(a4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(a5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(a6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(a7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
# ggplot(y) +
#   geom_hline(yintercept=0, lty="11", colour="grey30") +
#   geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
#   geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
#   geom_point(aes(cond, estimate)) +
#   labs(colour="", y = "Coefficient", x = "Condition") +
#   facet_grid(cols = vars(out)) + 
#   labs(title = "Coefficient Plot: Baseline OLS Model",
#        caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank()) 

y1 <- y
y1$mdl <- "Baseline"


### Demographic OLS Models #####################################################
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

## The Models
b1 <- lm(as.formula(paste("therm", iv1)), data = w)
b2 <- lm(as.formula(paste("p1.wash", iv1)), data = w)
b3 <- lm(as.formula(paste("p2.dist", iv1)), data = w)
b4 <- lm(as.formula(paste("p3.mask", iv1)), data = w)
b5 <- lm(as.formula(paste("p4.pray", iv1)), data = w)
b6 <- lm(as.formula(paste("behave", iv1)), data = w)
b7 <- lm(as.formula(paste("clicked", iv1)), data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(b1, conf.int = TRUE)
x2 <- tidy(b2, conf.int = TRUE)
x3 <- tidy(b3, conf.int = TRUE)
x4 <- tidy(b4, conf.int = TRUE)
x5 <- tidy(b5, conf.int = TRUE)
x6 <- tidy(b6, conf.int = TRUE)
x7 <- tidy(b7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(b1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(b2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(b3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(b4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(b5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(b6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(b7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
# ggplot(y) +
#   geom_hline(yintercept=0, lty="11", colour="grey30") +
#   geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
#   geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
#   geom_point(aes(cond, estimate)) +
#   labs(colour="", y = "Coefficient", x = "Condition") +
#   facet_grid(cols = vars(out)) + 
#   labs(title = "Coefficient Plot: OLS Model with Demographic Controls",
#        caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank()) 

y2 <- y
y2$mdl <- "Demographic Controls"


### OLS Models All Controls ####################################################
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

## The Models
c1 <- lm(as.formula(paste("therm", iv2)), data = w)
c2 <- lm(as.formula(paste("p1.wash", iv2)), data = w)
c3 <- lm(as.formula(paste("p2.dist", iv2)), data = w)
c4 <- lm(as.formula(paste("p3.mask", iv2)), data = w)
c5 <- lm(as.formula(paste("p4.pray", iv2)), data = w)
c6 <- lm(as.formula(paste("behave", iv2)), data = w)
c7 <- lm(as.formula(paste("clicked", iv2)), data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(c1, conf.int = TRUE)
x2 <- tidy(c2, conf.int = TRUE)
x3 <- tidy(c3, conf.int = TRUE)
x4 <- tidy(c4, conf.int = TRUE)
x5 <- tidy(c5, conf.int = TRUE)
x6 <- tidy(c6, conf.int = TRUE)
x7 <- tidy(c7, conf.int = TRUE)

x1$out <- "Thermometer"
x2$out <- "Policy: Hand Washing"
x3$out <- "Policy: Distancing"
x4$out <- "Policy: Mask Wearing"
x5$out <- "Policy: Prayer"
x6$out <- "Would Click"
x7$out <- "Clicked"

z1 <- tidy(c1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(c2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(c3, conf.int = TRUE, conf.level = .9)
z4 <- tidy(c4, conf.int = TRUE, conf.level = .9)
z5 <- tidy(c5, conf.int = TRUE, conf.level = .9)
z6 <- tidy(c6, conf.int = TRUE, conf.level = .9)
z7 <- tidy(c7, conf.int = TRUE, conf.level = .9)

z1$out <- "Thermometer"
z2$out <- "Policy: Hand Washing"
z3$out <- "Policy: Distancing"
z4$out <- "Policy: Mask Wearing"
z5$out <- "Policy: Prayer"
z6$out <- "Would Click"
z7$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)
z4 <- z4 %>% rename(low = conf.low, high = conf.high)
z5 <- z5 %>% rename(low = conf.low, high = conf.high)
z6 <- z6 %>% rename(low = conf.low, high = conf.high)
z7 <- z7 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3, z4, z5, z6, z7)
x <- rbind(x1, x2, x3, x4, x5, x6, x7)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Thermometer", 
                           "Policy: Hand Washing", 
                           "Policy: Distancing", 
                           "Policy: Mask Wearing", 
                           "Policy: Prayer", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
# ggplot(y) +
#   geom_hline(yintercept=0, lty="11", colour="grey30") +
#   geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
#   geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
#   geom_point(aes(cond, estimate)) +
#   labs(colour="", y = "Coefficient", x = "Condition") +
#   facet_grid(cols = vars(out)) + 
#   labs(title = "Coefficient Plot: OLS Model with Controls - Republican Sub-Sample",
#        caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank())

y3 <- y
y3$mdl <- "All Controls"


### Grid Plot of All Models ####################################################
rm(q, x, x1, x2, x3, x4, x5, x6, x7, y, z, z1, z2, z3, z4, z5, z6, z7)

y <- rbind(y1, y2, y3)

y$mdl <- factor(y$mdl, 
                levels = c("Baseline", 
                           "Demographic Controls", 
                           "All Controls"), 
                ordered = TRUE)

ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out), rows = vars(mdl)) + 
  labs(title = "Coefficient Plot: OLS Models - Republican Sub-Sample",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

## Cleaning Up
rm(list=setdiff(ls(), "working"))

```

## Synthetic Index{#syndex}

In this section, we collapse 5 of our outcome variables (the thermometer score and the 4 policy variables) into a single outcome variable using principal component analysis.  We then use the scores from the first component, which explains roughly 52% of the variation in the outcome variables, as a new dependent variable, and re-estimate our ATE specification.  Results are consistent with those presented in the main manuscript.

```{R princomp, fig.width = 10, fig.height = 8}
# Synthetic Index of Attitudinal Outcomes ######################################
## PCA and Scores ##############################################################
pc1 <- princomp(~ therm + p1.wash + p2.dist + p3.mask + p4.pray, 
                data = working, na.action = na.exclude, cor = TRUE)

summary(pc1)
working$pc.score <- pc1$scores[,1] 
```

```{R synthetic_index, fig.width = 10, fig.height = 8}
## New OLS with PCA as Outcome #################################################
iv1 <- " ~ treat + gender + age + education + party.f + covid.threat + covid.tested + covid.friends"
iv2 <- " ~ treat + gender + age + education + news.source + news.follow + covid.source + covid.follow + party.f + covid.threat + covid.tested + covid.friends"

w2 <- working %>% 
  filter(behave > 4)

### Baseline OLS Models ########################################################
## The Models
a1 <- lm(pc.score ~ treat, data = working)
a2 <- lm(behave ~ treat, data = working)
a3 <- lm(clicked ~ treat, data = w2)

## Create Dataframe for Coefficient Plot
x1 <- tidy(a1, conf.int = TRUE)
x2 <- tidy(a2, conf.int = TRUE)
x3 <- tidy(a3, conf.int = TRUE)

x1$out <- "Attitudinal Outcomes"
x2$out <- "Would Click"
x3$out <- "Clicked"

z1 <- tidy(a1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(a2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(a3, conf.int = TRUE, conf.level = .9)

z1$out <- "Attitudinal Outcomes"
z2$out <- "Would Click"
z3$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3)
x <- rbind(x1, x2, x3)

q <- z %>% full_join(x)

q$out <- factor(q$out, 
                levels = c("Attitudinal Outcomes", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
# ggplot(y) +
#   geom_hline(yintercept=0, lty="11", colour="grey30") +
#   geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
#   geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
#   geom_point(aes(cond, estimate)) +
#   labs(colour="", y = "Coefficient", x = "Condition") +
#   facet_grid(cols = vars(out)) + 
#   labs(title = "Coefficient Plot: Baseline OLS Model",
#        caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank()) +
#   scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

y1 <- y
y1$mdl <- "Baseline"

### Demographic OLS Models #####################################################
rm(q, x, x1, x2, x3, y, z, z1, z2, z3)

## The Models
b1 <- lm(as.formula(paste("pc.score", iv1)), data = working)
b2 <- lm(as.formula(paste("behave", iv1)), data = working)
b3 <- lm(as.formula(paste("clicked", iv1)), data = working)


## Create Dataframe for Coefficient Plot
x1 <- tidy(b1, conf.int = TRUE)
x2 <- tidy(b2, conf.int = TRUE)
x3 <- tidy(b3, conf.int = TRUE)

x1$out <- "Attitudinal Outcomes"
x2$out <- "Would Click"
x3$out <- "Clicked"

z1 <- tidy(b1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(b2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(b3, conf.int = TRUE, conf.level = .9)

z1$out <- "Attitudinal Outcomes"
z2$out <- "Would Click"
z3$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3)
x <- rbind(x1, x2, x3)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Attitudinal Outcomes", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
# ggplot(y) +
#   geom_hline(yintercept=0, lty="11", colour="grey30") +
#   geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
#   geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
#   geom_point(aes(cond, estimate)) +
#   labs(colour="", y = "Coefficient", x = "Condition") +
#   facet_grid(cols = vars(out)) + 
#   labs(title = "Coefficient Plot: OLS Model with Demographic Controls",
#        caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank()) +
#   scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

y2 <- y
y2$mdl <- "Demographic Controls"

### OLS Models All Controls ####################################################
rm(q, x, x1, x2, x3, y, z, z1, z2, z3)

## The Models
c1 <- lm(as.formula(paste("pc.score", iv2)), data = working)
c2 <- lm(as.formula(paste("behave", iv2)), data = working)
c3 <- lm(as.formula(paste("clicked", iv2)), data = working)

## Create Dataframe for Coefficient Plot
x1 <- tidy(c1, conf.int = TRUE)
x2 <- tidy(c2, conf.int = TRUE)
x3 <- tidy(c3, conf.int = TRUE)

x1$out <- "Attitudinal Outcomes"
x2$out <- "Would Click"
x3$out <- "Clicked"

z1 <- tidy(c1, conf.int = TRUE, conf.level = .9)
z2 <- tidy(c2, conf.int = TRUE, conf.level = .9)
z3 <- tidy(c3, conf.int = TRUE, conf.level = .9)

z1$out <- "Attitudinal Outcomes"
z2$out <- "Would Click"
z3$out <- "Clicked"

z1 <- z1 %>% rename(low = conf.low, high = conf.high)
z2 <- z2 %>% rename(low = conf.low, high = conf.high)
z3 <- z3 %>% rename(low = conf.low, high = conf.high)

## Combine Results in a DF to build Coefficient Plot
z <- rbind(z1, z2, z3)
x <- rbind(x1, x2, x3)

q <- z %>% full_join(x)

q <- q %>%
  filter(term == "(Intercept)" | term == "treatEndorsement" | term == "treatNorms")

q$out <- factor(q$out, 
                levels = c("Attitudinal Outcomes", 
                           "Would Click", 
                           "Clicked"), 
                ordered = TRUE)

q$cond <- dplyr::recode(q$term,
                        "(Intercept)" = "Intercept",
                        "treatEndorsement" = "Endorsement", 
                        "treatNorms" = "Norms")

## Filter out Intercept
y <- q %>%
  filter(cond != "Intercept")

## Coefficient Plot
# ggplot(y) +
#   geom_hline(yintercept=0, lty="11", colour="grey30") +
#   geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.2) +
#   geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
#   geom_point(aes(cond, estimate)) +
#   labs(colour="", y = "Coefficient", x = "Condition") +
#   facet_grid(cols = vars(out)) + 
#   labs(title = "Coefficient Plot: OLS Model with Controls",
#        caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
#   theme_bw() +
#   theme(panel.grid.major = element_blank(), 
#         panel.grid.minor = element_blank()) +
#   scale_y_continuous(breaks = c(-0.2, 0, 0.2, 0.4, 0.6, 0.8), limits = c(-0.2,0.8))

y3 <- y
y3$mdl <- "All Controls"

### Grid Plot of All Models ####################################################
rm(q, x, x1, x2, x3, y, z, z1, z2, z3)

y <- rbind(y1, y2, y3)

y$mdl <- factor(y$mdl, 
                levels = c("Baseline", 
                           "Demographic Controls", 
                           "All Controls"), 
                ordered = TRUE)

ggplot(y) +
  geom_hline(yintercept=0, lty="11", colour="grey30") +
  geom_errorbar(aes(cond, ymin=conf.low, ymax=conf.high), width=0.075) +
  geom_errorbar(aes(cond, ymin=low, ymax=high), lwd = 1.15, width=0) +
  geom_point(aes(cond, estimate)) +
  labs(colour="", y = "Coefficient", x = "Condition") +
  facet_grid(cols = vars(out), rows = vars(mdl)) + 
  labs(title = "Coefficient Plot: OLS Models",
       caption = "Shown with 90% and 95% confidence intervals. \n Note: Model for ''Clicked'' outcome uses subset data (only petition-eligible respondents).") +
  theme_bw() +
  theme(panel.grid.major = element_blank(), 
        panel.grid.minor = element_blank())

## Cleaning Up
rm(list=setdiff(ls(), "working"))
```

## Anchoring Effects

Below, we perform an exploratory analysis to asses potential reasons for the disparity in treatment effects between our thermometer score and the mask policy outcome.

We find that responses to the first two questions in the policy battery are significant predictors of responses to the third question (the masking policy question). Additionally, when we regress the masking policy question on our treatment and control for the answers to the first two battery questions, we find that the effect of our treatment is not statistically significant. We conclude that we have suggestive evidence that is consistent with the literature on conditional order effects, and that this type of anchoring may explain why we find a stronger effect of treatment on the thermometer question than on the battery question. 

```{R anchoring_effects}
pols3 <- lm(p3.mask ~ p1.wash, data = working)
# summary(pols3)

pols4 <- lm(p3.mask ~ p1.wash + p2.dist, data = working)
# summary(pols3)

pols5 <- lm(p3.mask ~ p1.wash + p2.dist + treat, data = working)
# summary(pols5)

xmodels <- list(pols3, pols4, pols5)

cm = c("(Intercept)" = "Intercept",
       "p1.wash" = "Policy: Hand Washing (Item 1)",
       "p2.dist" = "Policy: Distancing (Item 2)",
       "treatEndorsement" = "Endorsement",
       "treatNorms" = "Norms")

modelsummary(xmodels, 
             title = 'Mask Wearing (Item 3)',
             coef_map = cm,
             stars = TRUE, 
             gof_omit = 'IC|Log|Adj')
```


